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1. SUMMARY 


The objectives of this work are to continue the application of direct numerical simulation 
techniques to combustor flows and to use the simulations to develop a better understanding of 
the effects of turbulence on reactions. This work has involved three major tasks. First, we 
have performed simulations of chemical reactions without heat release in three dimensions to 
extend the results of our previous successful simulations. Second, we have performed high- 
resolution two and three-dimensional simulations of chemical reactions with heat release. 
These results have been compared with our previous calculations to determine the effects of 
heat release on reaction rate and flow field evolution. Finally, we have performed a series of 
simulations of three-dimensional turbulence decay for Dr. R. G. Deissler at NASA Lewis. We 
have designed, written, and debugged a code capable of performing simulations of three- 
dimensional reacting flows, and this code has been delivered to Mr. Russ Claus and Dr. 
Deissler at NASA Lewis to enhance the center’s computational capability. The work on this 
contract has been an important component of the basic combustion research being conducted 
at the NASA Lewis Research Center.! 1 ) 

The extension of the three-dimensional work from Part I of this report! 2 ) included some 
comparisons with laboratory data which showed good agreement in such quantities as average 
reactant concentration, rms fluctuating reactant concentration, rms fluctuating product 
concentration, and concentration correlations. Since there was no ad hoc modeling or 
adjustable parameters in these simulations, they have given strong evidence that this direct 
numerical simulation technique can accurately represent much of the important physics of this 
problem. 

This technique was next extended to the study of reacting flows with chemical heat release. 
This problem was solved using both the fully compressible equations as well as an 
approximate set of equations that is asymptotically valid for low Mach number flows. These 
latter equations have the computational advantage that high-frequency acoustic waves have 
been filtered out, allowing much larger time steps to be taken. The two-dimensional 
simulation results showed that the rate of chemical product formed, the thickness of the 
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mixing layer, and the amount of mass entrained into the layer all decrease with increasing 
rates of heat release. 


The subsequent three-dimensional simulations with heat release showed that, in agreement 
with previous laboratory experiments, the heat release is observed to lower the rate at which 
the mixing layer grows and to reduce the rate atwhicK chemical products are formed. The 


baroclinic torque and thermal expansion in the mixing layer y/ere found to produce changes in 
the flame vortex structure that act to produce more diffuse vortices than in the constant 
density case, resulting in lower rotation rates of fluid elements. Previously unexplained 
anomalies observed in the mean velocity profiles of reacting jets and mixing layers were 

J 

shown to result from vorticity generation by baroclinic torques. The density reductions also 
lowered the generation rates of turbulent kinetic energy and the turbulent shear stresses, 
resulting in less turbulent mixing of fluid elements. 


Calculations of the energy in the various wave number modes showed that the heat release has 
a stabilizing effect on the growth rates of individual modes. A linear stability analysis of a 
simplified model problem confirmed this, showing that low-density fluid in the mixing region 
will result in a shift of the frequency of the unstable modes to lower wave numbers. The 
growth rates of the unstable modes decrease, contributing to the slower growth of the mixing 
layer. 


Finally, an overall assessment of the status of the direct numerical simulation approach to the 
study of turbulent, reacting flows was made. While the strong constraints of spatial and 
temporal resolution of these methods will not be overcome in the near future, a great deal of 
useful information about the physics of such flows is available from these simulations, and the 
development of better subgrid-scale models will permit the implementation of the large eddy 
simulation technique to flows in more complex geometries. 
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2. SIMULATIONS OF A THREE-DIMENSIONAL, COLD. REACTING MIXING LAYER 


The results of the numerical simulations begun in the previous contract period were further 
analyzed and compared with experimental data. The simulations agree very well with the 
predictions of self-similarity theory, yielding approximately linear growth rates of various 
computed length scales, including the mean velocity half-width, the mean vorticity thickness, 
and the mean product thickness. In addition, mean profiles of the average reactant 
concentrations, the rms fluctuating reactant concentrations, the concentration correlations, the 
average product concentrations, and the rms fluctuating product concentrations collapsed in a 
manner consistent with self-similarity theory. This was not an artifact of the initial condition?, 
since the initial concentration fields were not in the self-similar form, and the self-similarity 
was maintained over a time in which the width of the layer grew by a factor of 5. Simulation 
results also compared reasonably with laboratory data. Computed profiles that were 
qualitatively similar to corresponding laboratory profiles were obtained for the average 
reactant concentrations, the rms fluctuating product concentrations, and the concentration 
correlations. Computed profiles of average product concentrations were in approximate 
agreement with laboratory data. 

In the course of this work, a detailed analysis was made of the accuracy of these techniques 
when applied to the simulation of chemically reacting flows. When the reaction rate is large 
relative to the species difTusivity, steep gradients can develop in the concentration fields 
during the course of the simulation, even though the initial fields are smooth. Thus, the 
accuracy of a simulation can decrease with time until significant errors develop. Some test 
calculations were performed in which the exact analytical solution was known so that the 
resulting numerical errors could be quantitatively estimated. These simulations showed that in 
the limit as the time stepping errors approached zero, the superalgebraic convergence (or 
infinite order accuracy) characteristic of spectral methods could be maintained. 

The results of this work are described in much greater detail in the two papers reproduced in 
two papers by Riley et al.! 3 I M 
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3. TWO-DIMENSIONAL MIXING LAYERS WITH HEAT RELEASE 

The purpose of this aspect of the work was to begin an investigation into the interactions 
between chemical heat release and fluid dynamics in a mixing layer. This is an extension of 
the previous work in which the chemical reaction was a passive process, and thus did not 
affect the fluid motion. In a reacting flow with heat release, the dynamics of the fluid motion 
are coupled with the chemical reaction through the inhomogeneous density distribution 
caused by the thermal expansion. This problem can be attacked by solving the compressible 
Navier-Stokcs equations with heat production, together with the species transport equations. 
This set of equations contains the vorticity, entropy and acoustic modes, which can be of 
greatly different frequencies. In particular, for low subsonic fluid motion, the acoustic modes 
are in frequency bands much higher than the other two modes. Due to this high frequency, 
the acoustic fluctuations do not interact effectively with either the vorticity mode or the 
entropy mode. From the computational point of view, tracking the acoustic fluctuations 
requires extremely small time steps, thus decreasing the efficiency of the computation. To 
mitigate this constraint, a set of approximate equations was derived, asymptotically valid for 
small Mach number flows. 

The validity of the low Mach number approximation in the parameter range of the simulations 
performed here was tested by solving the exact, fully compressible equations for a reacting 
mixing layer with a Mach number of 0.2. Vorticity contours obtained from solutions of the 
exact and approximate equations showed that although acoustic waves propagate throughout 
the flow field, they apparently have no influence on the vorticity dynamics. Simulations of the 
two-dimensional mixing layer undergoing vortex rollup showed that the effect of heat release 
was to reduce significantly the amount of product produced in the layer. Because of the 
decrease in density caused by exothermic chemical reactions, the beat release and resulting 
density changes act to decrease the rate of mass entrainment into the layer, since the overall 
layer growth is approximately unchanged from the cold flow case. This is in qualitative 
agreement with the experimental findings of Wallace,^ who conducted an experimental study 
of mixing layers that included a chemical reaction with weak heat release. 

The variations in density permit a major new source of vorticity generation in the form of 
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baroclinic torques. This tends to reduce the vorticity near the outer edges of the layer, 
inhibiting the rollup process. Another effect of heat release is to raise the viscosity of the 
flow. This can also have a stabilizing effect on the flow. 

An important result of this aspect of the work was to demonstrate the effectiveness of the 
direct numerical simulation approach as a predictive tool in studying complex reacting flows 
with heat release. This clearly suggested the importance of extending this computational 
capability to three-dimensional flows. A detailed discussion of the derivation of the low Mach 
number approximation equations, the numerical solution procedure and results of the 
simulations are given in the paper reproduced in Appendix A.I 6 1 
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4. HEAT RELEASE IN CHEMICALLY REACTING TURBULENT MIXING LAYERS 

The results of the two-dimensional simulations with heat release showed that the presence of 
heat release produced a stabilizing effect on the flow field. However, in fully turbulent mixing 
layers, other secondary modes of instability are present which could change this situation. In 
order to investigate this possibility, a fully three-dimensional code was developed to simulated 
reacting flows under the following conditions. First, to ease time stepping stability constraints, 
low Mach number approximation equations were used to filter out acoustic waves. Second, a 
single binary, irreversible reaction between two species to form a product was used. The 
reaction rate was independent of the temperature. Third, the viscosity, thermal and molecular 
diffusivities and the specific heats were taken to be temperature- independent constants. 
Finally, the amount of heat release was restricted in order to ensure that large velocities would 
not be generated by the reaction to violate the low Mach number approximation. 

The computational domain was chosen to be large enough to contain the most unstable mode 
and its subharmonic. The velocity field was initialized by adding a hyperbolic tangent 
velocity profile to a low- level, three-dimensional, broad-band background perturbation field. 
In addition, in the "forced" simulations, an additional perturbation in the form of the most 
unstable mode and its subharmonic was added. These forced simulations are analogous to 
laboratory experiments in which well-defined harmonic perturbations are introduced into the 
flow at, or upstream of, the splitter plate. The chemical reactant fields were initially two- 
dimensional and consisted of two species fields that were offset so that there was initially no 
overlap between them. 

In three spatial dimensions, an additional secondary instability mechanism comes into play. 
At large amplitudes, it is manifest as streamwise, counter-rotating vortex tubes on the braids 
between the vortex cores. The effect of these structures on the flow field is to induce a 
velocity field that acts to pump fluid between the two layers, thereby further convoluting the 
reaction interface. This tends to increase mixing between the two streams and enhance the 
chemical reaction rate. 

In spite of these additional instabilities, the effect of heat release on the chemical product 
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formation is very much similar to the two-dimensional results. Heat release tends to reduce 
the product generation, again by inhibiting the mixing between the two streams. When the 
flow is forced, by including the two-dimensional most unstable mode and its subharmonic, the 
product generation is greater than without forcing but is still significantly less than in the 
comparable run without heat release. However, the three-dimensional modes do enhance the 
product generation compared to the strictly two-dimensional runs by providing an additional 
mechanism for wrinkling the flame front and increasing the intercomponent mixing. For the 
growth of the mixing layer thickness, there is a small initial enhancement, followed by a 
decrease in the layer growth rate. These results are qualitatively similar to those obtained 
experimentally by Wallace^ 5 ! and HermansonW and are consequences of lower rates of flujej 
entrainment into the mixing region when exothermic chemical reactions occur. 

An analysis of the turbulent kinetic energy equation for these simulations showed that the 
most important effects of heat release are significant reductions in the turbulence production 
and turbulent transport. In the case of the production term, the lower density resulting from 
the combustion is directly responsible for most of this change. Transport by the fluctuating 
motion decreases in the heat release runs simply because the turbulence levels are lower. The 
heat release results in a production of turbulent kinetic energy along the center of the mixing 
layer through the expansion part of the velocity-pressure gradient correlation. However, the 
production is small in these simulations compared with the decrease in the exchange of energy 
with the mean flow, yielding an overall lower turbulent kinetic energy profile. As heat release 
is increased, there can be a significantly greater conversion of internal energy to kinetic 
energy. With much higher heat release, this could possibly result in an overall increase in the 
kinetic energy. 

The reduction in mixing layer growth and turbulence production observed in the simulations 
with heat release is also consistent with the results of a simplified linear analysis based on 
linear profile approximations to the mean velocity density fields. This analysis shows that, iq 
the presence of heat release, there is a shift of the most unstable wave number to a lower value 
and a decrease in the value of the largest wave number that is still unstable, and there is also a 
decrease in the growth rate of the most unstable mode. The precise effect of these stabilizing 
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factors depends on the relative time scales involved in the manifestation of the primary 
instabilities and the time scales over which the density decreases develop. However, the net 
result of this analysis is to show that heat release can be strongly stabilizing, both in the linear 
as well as in the nonlinear regime. 

A more complete analysis of these simulations is given in a paper by McMurtry et alJ 8 l and the 
paper provided in Appendix B,l°l as well as in the Ph.D. thesis of McMurtryJ 10 ! 
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5. FURTHER ANALYSIS OF THREE-DIMENSIONAL SIMULATIONS 
5.1 Computation of Probability Density Functions 

We have employed results of direct numerical simulations to construct the probability density 
function of a conserved scalar variable in a perturbed three-dimensional temporally evolving 
mixing layer. Calculations have been performed to a time just before mixing transition, and 
the results are comparable to experimental data obtained in the same region. Some additional 
simulations were performed to explore the possibility of extending the simulation procedure to 
larger domains. 

Probability density functions have proven useful in the theoretical treatment of turbulent 
reacting flows since the early work of Hawthorn et alJ 11 ) An approach based on the single- 
point probability density function of the scalar variables has the advantage that the effects of 
chemical reactions appear in closed form, eliminating the need for any turbulence modeling 
associated with the scalar-scalar correlations. However, models are needed for the closure of 
the molecular mixing term and also the turbulent convection. 

In previous work,I 12 l l 13 l I 14 l the pdf approach was employed in a variety of turbulent reacting 
shear flows, and the results were compared with available experimental data under similar 
hydrochemical conditions. In this work, the effects of the molecular mixing term were 
modeled by a coalescence/dispersion model, and the turbulent flux of the pdf was modeled by 
a simple gradient diffusion approximation. 

As far as the first few moments are concerned, the results obtained from the pdf transport 
equation are in reasonable agreement with experimental data. However, a major difference 
between the calculated and measured shape of the pdf was observed in the two-stream plane 
mixing layer calculations. 

The experimentally observed shape of the pdf is dependent on the streamwise location, where 
the pdf is measured. Measurements of Konradl 16 ! and Koochesfahanil 16 ! in nonreacting mixing 
layers indicate that below the mixing transition point, the concentration field consists of pure 
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unmixed fluid originating from either the high- or low-speed stream, and mixed fluid is found 
only at the thin interface separating the two fluids. Under these conditions, the shape of the 
pdf is rather simple. It has spikes only at values corresponding to the free stream 
concentrations and is negligible at other concentration values. Above the point of mixing 
transition, the experimental measurements indicate that there is a change in the shape of the 
pdf. It still consists of the two spikes at the concentration values corresponding to the free 
stream concentrations. However, a third peak in the mixing region of the shear layer is 
observed (Figure 1). The asymmetry of this peak indicates the excess of high-speed fluid 
species in the mixing core of the layer. Furthermore, experimental measurements show that 
the functional form of the pdf (including the position of the third peak) changes very little as 
the mixing layer is transversed. 

Calculations based on a modeled pdf transport equation by Givi et alJ 14 l do not predict the 
shape of the pdf consistent with experimental observations. The predicted results indicate 
that the location of the pdf peak, with respect to the concentration coordinate, changes as the 
layer is transversed. The major reason for this discrepancy is due to the shortcomings 
associated with the gradient diffusion modeling of the turbulent flux of the pdf. It has been 
mentioned by Givi et alJ 14 l that in a mixing layer the inner turbulent region is often 
interrupted by the presence of irrotational surrounding flow. Therefore, a simple gradient 
diffusion assumption, which is a reasonable model for fully turbulent flows, would not be 
expected to result in plausible predictions in a highly intermittent flow such as the mixing 
layer. 

With direct numerical simulation, it is possible to simulate the mixing layer directly without 
any need for additional modeling, as is usually required if the equations are statistically 
averaged. The results of the direct numerical simulations can be analyzed in order to obtain 
useful statistical information about the important vector and scalar variables. Comparison of 
such statistical quantities with experimental data is very useful in obtaining a better 
understanding of the important physical phenomena that occur in turbulent flows. The 
probability density function of a conserved scalar is constructed directly from the flow fields 
generated by the direct numerical simulations. The following results have been obtained for 
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the prc-transitional region of the mixing layer. 


The initialization of the mean flow and the velocity perturbation fields is consistent with the 
RUN R simulations reported by Metcalfe et alJ 17 ! The initial small amplitude random 
perturbations are superimposed on a mean hyperbolic tangent velocity profile. These linear 
perturbations grow as time progresses, and, as a result, the initially smooth flow experiences a 
period of transition before it reaches a fully turbulent state. 


One way of determining whether the flow is in the pre- or post-transitional domain is to 
examine the growth of the vorticity thickness and the product thickness. Vorticity thickness 
is defined by 


Su) = 


A U 

,dU) 

\ /max 


The product thickness is a measure of the degree of mixing in the core of the mixing layer. 
One way of measuring the degree of mixing, experimentally, is to inject reactants on two sides 
of the mixing layer and measure the amount of product resulting from the chemical reaction 
between the two reactants. The reactants selected must have very fast chemical kinetics so 
that the amount of product formed is proportional to the amount of mixing. In such a system, 
the product thickness is defined by 


6p - 


( C B) a 


f Cpdy 


where C p is the concentration of product and (Cb) m refers to the concentration in one of the 
streams. 


Experimental evidence! 18 ! indicates that, before mixing transition, both product thickness and 
vorticity thickness increase in the streamwise direction, while the ratio between the two 
remains more or less constant until the transitional point. Beyond this point, there is an 
abrupt increase in product formation and, therefore, a sharp increase in the ratio of the 
product thickness to the vorticity thickness. In one simulation, the product thickness and the 
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vorticity thickness, as well as the ratio between the two, are calculated directly from the 
instantaneous velocity and scalar fields. The product concentration is calculated by 
employing a fast chemistry model (flame sheet approximation! 13 !). This approximation relates 
the product concentration to the concentration of a conserved Shvab-Zeldovich scalar 
variable. 

A plot of the normalized product thickness, normalized vorticity thickness, and the ratio 
between the two is presented as a function of time in Figure 2. Note that both thicknesses 
increase while the ratio between the two remains constant. At the time of t * = 160, boundary 
effects become important, so this calculation was not continued beyond this point. The flow 
at this time is still in the pre-transitional mixing region, since there has not been a sharp 
increase in the ratio of the product thickness to the vorticity thickness. 

The probability density function of a conserved scalar variable at the time of t * = 1 60 is 
presented in Figure 3. The instantaneous free stream value of the concentration is equal to 
zero in the top stream and equal to 1 in the bottom stream. Analysis of this figure suggests 
that the mechanism of mixing before mixing transition is fairly simple in that the fluid across 
the shear layer width is essentially composed only of the fluid originating from the two 
streams. The experimental measurements of Koochesfahani! 16 ! also show the same behavior 
for the pdf in the pre-transitional region. In this region, the fluid is either from the top stream 
or from the bottom stream, and strong mixing has not yet occurred. 

5.2 Simulations on Larger Domains 

In an attempt to carry the above simulations beyond the mixing transition regime, some 
studies were made as to the feasibility of reinitializing the code at the pre-transition point but 
on a larger domain. The basic procedure involved doubling the size of the computational 
domain in all three dimensions. The regions above and below the mixing layer were extended 
uniformly based on the values of the computational variables at the top and bottom 
boundaries. In the streamwise and spanwise directions, the periodicity length was doubled 
and the initial data were prescribed by retaining only alternate data points in the original 
calculation but repeating the sequence twice. Both calculations were then continued in time 
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so that they could be compared. 


Figure 4 shows a three-dimensional perspective plot of the species concentration field C B at a 
time of / * = 90. This is after the first pairing and at the beginning of the second pairing on a 
domain of side 8jt. The state of this field at / * = 120 is shown in Figure 5, at which time the 
second pairing is taking place. The corresponding plots for the larger 16rr domain are shown 
in Figures 6 and 7. This computation is then carried on further to t* — 150 as shown in 
Figure 8. 

While the perspective plots indicate strong similarities between the original and expanded 
simulations, a more quantitative analysis can be undertaken by examining the evolution of 
some of the low-order modes in the two simulations. Figure 9 shows the evolution of the 
energies in the mean flow and the second and third subharmonics for the original simulation. 
Note that the second subharmonic {E sa 3 ) saturates by about /*=80, while the third 
subharmonic saturates by about t * = 130. The behaviors of these modes in the simulations on 
the expanded domain from / * = 90 to / * = 150 are shown in Figure 10. Note that while there 
are some differences, the basic trends are similar, especially considering the fact that 
boundary effects are becoming significant in the original run by / * = 150. 

Based on these simulations, this grid expansion technique appears to have significant potential 
as a tool to further the study of reacting mixing layer dynamics well into the mixing transition 
regime. 
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6. CONCLUSIONS AND DISCUSSION 


This contract research work has provided the first major applications of the direct numerical 
simulation approach to reacting flows. In the first place, the possibility of employing high- 
resolution spectral numerical methods to reacting flows was examined and shown to be 
feasible through a series of test calculations. That the retention of the critical superalgebraic 
convergence property of spectral methods was possible when applied to the simulation of 
moderately fast reactions in a mixing layer undergoing nonlinear vortex pairing was 
demonstrated quantitatively. The three-dimensional simulations showed the relative 
importance of secondary instabilities in the flow and the role they have in mixing 
enhancement. Comparisons of the results of these simulations of a fully turbulent flow with 
experiments and self- similarity theory showed excellent agreement on a number of statistical 
quantities, especially considering that there were no adjustable parameters or ad hoc constants 
in the model. 

The extension of this work to the more complex case with heat release was successfully 
undertaken. Initial simulations, involving the pairing of large-scale, two-dimensional vortical 
structures in the mixing layer, showed that heat release can have a strongly stabilizing effect 
on the flow, in agreement with laboratory experiments. An analysis of these and later three- 
dimensional simulations showed that this stabilization was due to both linear and nonlinear 
effects. The thermal expansion tends to reduce the growth rates of the unstable modes in a 
laminar mixing layer, as does the reduction in Reynolds stress production and generation of 
baroclinic torques in the nonlinear regime. 

The principal advantages of the direct numerical simulation approach include the following: 
the flow can be examined in detail, since all quantities are known at each point in space and 
time; parameters can be easily varied and experimental conditions are easily controlled; large- 
scale structures are directly computed. The main disadvantage is in the limited spatial and 
temporal resolution available. For reacting turbulent flows at realistic Reynolds and 
Damkohler numbers, it is often necessary to include a subgrid-scale model. The development 
of accurate subgrid-scale models for such simulations [often referred to as large eddy 
simulations (LES) to distinguish them from the unmodeled full turbulence simulations (FTS)] 
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is one of the most critical issues that needs to be addressed in order to extend the power of 
these methods to more practical situations. Present subgrid-scale models for nonreacting flows 
rely on hypotheses regarding the nonlinear cascade of energy from large to small scales and 
the universality of motion at high wave numbers. For large Damkohler numbers, heat release 
and other chemical reaction effects could significantly modify such behavior. Thus, the 
development of accurate, physically derived subgrid-scale models is a very challenging task. 
One particularly promising approach is the application of renormalization group methods.! 10 ! A 
more detailed comparison of the various approaches to the simulation of reacting turbulent 
flows is given in a paper by Jou and Riley. I 20 ! 

The simulations performed in this research work indicate that while there are limitations on 
Reynolds number and Damkohler number, important aspects of the physics are being properly 
represented and that the velocity and concentration fields being generated by these 
simulations provide a useful source of information about the dynamics of the reacting flow. 
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Figure 1. Comparison of measured and calculated pdfs of a normalized conserved scalar (the 
mixture fraction, 0 at points across the mixing layer. The similarity coordinate, is the ratio 
of the cross-stream coordinate divided by the downstream distance from the virtual origin of 
the mixing layer. These pdfs represent conditions at a distance of « 15 cm from the virtual 
origin. 



18 









Figure 4. Three-dimensional perspective plot of species concentration field C B at t * = 90 on 
8?r domain. 
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Figure 5. Three-dimensional perspective 


plot of species concentration field C B at t - 120 on 


8jt domain. 
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Figure 6. Three-dimensional perspective plot of species concentration field C B at t ' = 90 on 
16 tt domain. 
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Figure 7. Three-dimensional perspective plot of species concentration field C* at = 120 on 
16* domain. 
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Figure 8. Three-dimensional perspective plot of 
1 6 rr domain. 


species concentration field C B at l* = 150 on 
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Figure 9. The evolution of the energies in the mean flow, E M , and the second and third 
subharmonics, E S h 2 ,E S h 8 , for the turbulent mixing layer simulation in the 8tt domain. 
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Figure 10. The evolution of the energies in the mean flow, Em, and the second and third 
subharmonics, E S H a ,£su a , for the turbulent mixing layer simulation in the \6n domain. 
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APPENDIX A 


DIRECT NUMERICAL SIMULATIONS OF A REACTING 
MIXING LAYER WITH CHEMICAL HEAT RELEASE 


P.A. McMurtry*. W.-H. Jou,** 
J.J. Riley'*’, and R.W. Metcalfe+t 
Flow Industries, Inc. 

Kent, Washington 


Abstract 

In order to study the coupling between chemical, 
best reltsse and fluid dynasties, we have perforned 
direct numerical simulations of a chemically react- 
ing mixing layer with heat release. We have treated 
the fully compressible equations es well as an ap- 
proximate set of aquations that is ssymptotically 
valid for low Mach number flows. These latter 
equations have the computational advantage that high 
frequency acoustic waves have been filtered out, 
allowing much larger time steps to be taken in the 
numerical solution procedure. A detailed derivation 
of these equations along with an outline of the 
numerical solution technique is given. Simulation 
results indicate that the rate of chemical product 
formed, the thickness of the mixing layer, and the 
amount of mats entrained into the layer all decrease 
with increasing rates of heat release. 


nomenclature 

Cj - molar concentration of ith chemical species 
Ce - nondimens ional heat release parameter 
C v - apecific heat at constant volume 
Da - Damkolher number 
e - internal energy 

1 0 - wavelength of fundamental perturbation a»de 
M - Mach number 
p - pressure 

p° - thermodynamic pressure 
pi - dynamic pressure 
P - dimensionless heat source 
Pe - Peclet number 
- heat flux vector 
Q - rate of heat release 

Rj - rate of reaction of ith chemical species 
Re - Reynolds number 
t - time 
T - temperature 

U - velocity difference across mixing layer 
v - magnitude of velocity 

^ - velocity vector 

$ - position vector 

Y " ratio of specific heats 

c - ym 2 
p - density 
T - stress tensor 


1. Introduction 

The development of accurate combustion models 
relies on understanding more fully the fundamental 
interactions among the thermal, chemical, end fluid 
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dynamical processes in reacting flow systems. The 
work reported here is directed at studying the 
interactions between chemical heat release and 
fluid dynamics, one of the least understood aspects 
of combustion. To achieve this purpose we have 
performed direct numerical simulations of a two- 
dimensional, temporally growing mixing layer, 
including a single step, irreversible chemical 
reaction between initially unmixed chemical species. 
The idea of direct numerical simulations is to use 
very accurate and efficient numerical methods to 
solve the governing equations for the detailed time 
evolution of the complex flow field. The appeal of 
this method is that no closure modelling is neces- 
sary and the complete flow field is known at every 
instant ao that any statiatical quantitiea of 
interest may be computed. The main disadvantage of 
this method is the limited range of spatial and 
temporal scalea that can be resolved, restricting 
accurate numerical resolution to moderate Reynolds 
numbera. No turbulence or transport modelling has 
been attempted, limiting our study to the large 
scale features of the fluid motion. However, the 
large-scale motions play a significant role in the 
evolution of a mixing layer and do not depend 
Strongly on the Reynolds number 1 . 

Previous direct numerical simulations of an 
incompressible reacting mixing layer without heat 
release 2 focused on how the turbulent flow field 
affects chemical reaction rates. This work has 
given physical insight into how the vortex rollup 
and three-dimensional turbulence affect the chemi- 
cal reaction. Furthermore, from comparison of 
simulation results with laboratory data, this work 
hat also led to increased confidence in the appli- 
cation of the direct numerical simulation methodo- 
logy to this class of problems. However, because 
no heat was released, the chemictl reaction was a 
passive process and did not influence the fluid 
motion. In a reacting flow with heat release, the 
dynamics of the fluid motion are coupled with the 
chemical reaction through the nonhomogeneous den- 
sity distribution esused by the thermal expansion. 
The work presented here attempts to include coupling 
between the chemical reaction, the heat release, end 
the flow field. 

To accomplish this goal, one may proceed to 
solve the compressible Kavier-Stokes equation with 
heat production, together with the species transport 
aquations. This set of aquations contain the vorti- 
city, entropy and acoustic modes from the linearised 
disturbance point of view.3 In many practical 
cases, these nodes are of greatly different fre- 
quency. Particularly, for lov subsonic fluid 
motion, the acoustic modes are in frequency bands 
much higher than the other two i»des. Since the 
frequency of these modes is high, the acoustic 
fluctuations do not interact effectively with either 
the vorticity mode or the entropy node. From the 
computetional point of view, tracking the acoustic 
fluctations requires extremely email time steps, 
thus decreasing the efficiency of the computation 
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as far as the vortex/entropy dynamic* are concerned. Governing Equations 
For these reasons, we have derived a aet of ap- 
proximate equations asymptotically valid for small The conservation equations of mats, momentum, 

Kach cumber flows. This derivation is given in energy, and chemical species may be written in the 

Section 2. following form. 


In Section 3, a numerical solution procedure to 
solve the governing equations is outlined. The 
differences between the method used here and 
similar solution techniques for unsteady, imcom- 
prcssible fluid flows are pointed out. 


Although a turbulent mixing layer ia inherently 
three-dimensional, the two-dimensional simulations 
discussed here reveal important information about 
eome of the dominant features. Laboratory experi- 
ments have shown that, at least in its early stages, 
the mixing layer ia dominated by large-scale, two- 
dimensional vortices, with the growth of the lsyer 
dominated by the pairing of these vortices.^* Two- 
dimensional simulations have been shown to accur- 
ately portray the characteristic large-scale rollup 
and vortex pairing process of mixing layers. 5 
Implications of heat releasing chemistry on this 
process can thus be meaningfully addressed with 
two-distensional simulations. In Section A, simula- 
tion results with and without heat releasing 
chemistry are presented. 


■|£- ♦ V*(pv) ■ 0 (la) 

pg--V P fH (lb) 

p E ■ - V*(pv) - V*q ♦ V*(T*v) ♦ Q (lc) 

3C. 

♦ V'f^v) - R £ 4 Va.VCj) , (Id) 

where 


In a perfect gas in which the molecular weights of 
the individual species sre approximately equal, the 
equation of atate is 

p ■ pRT . (le) 


2. Low Mach Humber Approximation 

Detailed numerical solutions to general combus- 
tion processes are prohibited by exceasive computer 
storage requirements and unacceptable computational 
expense. This results in part from the wide range 
of time and space scales present in reacting flow 
problems. For example, accurate resolution of steep 
grsdients of reacting species msy require a grid 
spacing orders of Mgnitude finer than that neces- 
sary to resolve other characteristic flow struc- 
tures. Similarly, time scales characterising acous- 
tic wave propagation may be orders of magnitude 
smaller than time acales characterising convection 
processes. It is this last problem to which atten- 
tion ia directed here. 

The existence of high frequency acoustic waves 
places a severe restriction on the time stepping 
increments used to numerically advance the fully 
compressible equations in time. To deal with this 
difficulty an approxinate set of equations has been 
derived for low Kach number flows. These aquations 
have the computational advantage that acoustic 
waves have been filtered out, relieving the above 
mentioned time stepping constraint. However, at 
the same time they allow density nonuniformities 
resulting from heat rcleaae to develop. In the 
following, we have atarted from the exact governing 
equations. By using a small Mach number expansion, 
a aet of ordered approximate equationa ia obtained. 

A set of criteria on other non-dia*nsional para- 
meters for the problem can also be obtained for the 
applicability of the approxisution. The approach 
used is similar to that of Rehm and Baum,^ who 
derived approximate equations for buoyantly driven 
flows with slow heat addition. Application of 
similar equations in subsonic reacting flows have 
been discussed by Oran and Boris^, although the 
distinction between the lowest order pressure and 
the first order pressure was not pointed out. This 
distinction ia essential in both a theoretical 
framework and the numerical procedure. 


The variables are scaled by their respective refer- 
ence quantities as 

P - P 0 P' 
v - D v* 


P 

e 

Q 

«* 

X 


t 


(PoRToV 


< C vV e ’ 


V 

L x' 
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( 2 ) 


In (2), the primed variables are dimensionless and 
the subscripted variables refer to the reference 
state. The resulting nondimensional equations, 
with the primes now dropped from the dimensionless 
variables, are 


If* V*(f£) - 0 

(3a) 

v* 2 p " - y p ♦ £r v* 

(3b) 

P jL E . - VU^DpS) - ~|L 

♦ < v *(?*v)) ♦ CeP 

(3c) 

8C. 

7T* V * (C i^ "Da R. vijV 2 C. 

(3d) 

P • PT 

(3e) 

e - t ♦ Y(-r-i) mV . 
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To close this system of equations the folloving 
constitutive relations Bust be provided 


R i “ R i (C i* T) 
P - POO 

q " q(VT) 
t ■ T(Vv) . 


. Pr, Re, Da, and Pe are the Prandtl, Reynolda, 
Damkb'hler, and Peclet numbera, respectively. The 
parameter Ce is given by 


Ce * 


Q n L 

o o 

D p C T 
o o V o 


For the validity of the low Mach number approxima- 
tion, Q 0 Bust be restricted such that Ce ■ o(l). 

If Ce » 1, large velocities will be generated by 
the combustion, and the low Mach number approxima- 
tion will not be valid. 


Small Mach Number Expansion 


Defining c ■ tM^ and assuming that c «1, the 
dependent variables may be expanded as power series 
in t: 


„ . JO) „ r Jl) . 

p ■ p ♦cp ♦ . . . 

P - p (0) 4 tpd) 4 ... 

T - T (0) ♦ CT (1) 4 ... 
C i “ C [ 0) * EC [ 1> 4 “• 


(4) 


Substituting equations (4) in (3) and collecting 
the teroth order terms in c results in the 
folloving lowest order equations: 


.jo) r 
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T> 

O 

o 

> 

o 

(5a) 

Vp (0) - 0 


(5b) 

JO) D (0) 

p Dt - T 

(0) . . (rl)p C0) ? 4(0) 

4 db v ^ (0) 4 CeR 

(5c) 

8c( i 0) r 

c| 0) : (0) l - Da R| 4 f V 2 C| 0) 
i J i Pe i 

(5d) 

,<»> . „<•> 

T (o) . 

(5e) 


It it observed from this last set of equations that 
the momentum equation haa degenerated to describe 
the variation of the thermodynamic preasure p*0). 
To complete the description of the velocity field, 
$'0) ( the first order momentum equation must be 
included and is given as 


(0) ^. 0) *W m 

* Dt 




$0) 


(5f ) 


It should be noted that all dependent variables 
appear only to lowest order except pressure, in 
which both p'°) an d p'D appear. The superscripts 
will now be dropped except to distinguish p'®) and 
pd'. The lowest order momentum equation eimply 
states that in the first approximation the thermo- 


dynamic pressure, p^®), is constant in space but 
auy vary with tism. In this approximation the 
sound speed it infinite ao that any disturbances in 
thermodynamic preasure caused by combustion are 
felt instantaneously throughout the fluid. For 
combustion in an open domain, p^®) * po* and the 
combustion is a constant preasure process. p'D is 
the dynamic pressure associated vith the fluid flow 
and, to the lowest order, does not participate in 
thermodynamic processes. 


These equations may be stated in an alternate 
form that proves to be helpful in obtaining numeri- 
cal solutions. Equations (5a) and (Sc) can be 
combined to give 


V*v 




PrRe dt 


* CeP 


(I 


( 6 ) 


In a closed domain vith adiabatic walls or in a 
domain with periodic boundary conditions, equation 
(6) can be integrated to give 

Under these conditions the combustion is a constant 
volume process and the pressure in a closed domain 
increases due to the rate of chemical energy re- 
leased in the chamber. Equations (6) and (7) are 
used to replace equations (5b) and (5c), giving the 
final set of equations. Note that if the heat 
release term, t “ 0, then V*$ ■ 0 and these equa- 
tions reduce to those of a strictly. incompressible 
fluid. 


3. Numerical Solution Procedure 


Equations (5a) and (5f) may be written in the 
following form: 

9 


* pV*v 4 v*Vp ■ 0 
^ (p$) - - Vp (1) ♦ A(x) 


( 8 ) 

(9) 


Here A(x) represents the viscous and convective 
components of the momentum equation. The local 
density is obtained by aubetituting equation (6) 
into (8) and advancing the density ahead in time 
using a second order Adams-Bashforth time stepping 
scheme. The chemical species equations are also 
solved using an Adams-Bashforth scheme. All 
spatial dcrivitives are computed using pseudo- 
spectral methods. The preasure and velocity 
components are obtained by eolving equation (9) in 
two stages. First equation (9) is advanced to an 
intermediate time step neglecting preasure effects: 

pv* - pv 5 " 1 4 AtJj A(x) n - ^ A(x)* 1-1 ] . (10) 

The remaining preasure term gives: 

Vp ci )" n/2 (11) 

Note that adding equations (10) and (11) results in 
a second order accurate finite difference approxima- 
tion to the time derivative of P$. Taking the 
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divergence of equation (11) give* the following 
Poieson equation for preaaure: 


„ 2 (ir 

-V p 


V*(pv) T 


V‘(pv)* 


Since we are treating a spatially periodic problem 
in thie work (aee Section 4), equation (12) ia 
eaaily aolved for p' 1 ' by taking the Fourier trans- 
form of aquation (12), aolving for the tranaform of 
p^l), and then transforming back to physical apace. 
The velocity cooponenta are obtained by uaing thia 
value of p'*' in equation (11). This sKthod ia sim- 
ilar to nuaierical solution techniques for the incom- 
pressible, unsteady, Navier-Stokes equations as out- 
lined by Peyret and Taylor®. The important differ- 
ence is that the density is not constant and an ap- 
proximation for V*(p$) n ^^ must be obtained directly 
from equation (5a) using computed values of p** 4 *, 



4. Simulation Results 

In this section we discuss some preliminary 
results of our numerical simulations of the temporal 
development of a chemically reacting mixing layer 
with heat release. The computational domain ia a 
rectangle with 64 x 64 equally spaced grid points. 
The flow is assumed periodic in the streamwise 
direction and free-slip, impermeable boundary condi- 
tions are employed at the transverse boundary. The 
initialisation of the flow field ia the same as 
used by Riley and Metcalfe^ in their two-dimensional 
simulations of an incompressible reacting mixing 
layer without heat release. Namely, the velocity 
field was initialised with a hyperbolic tangent pro- 
file, and perturbations corresponding to the most 
unstable mode and its subharmonic as determined from 
linear stability theory^. The chemical reactant 
concentrations were set at follows: 

y 

c i ( ** o) ‘ 777r / «p<-c 2 /yo><ic 

0 Zee 

C^tx.o) - 1 - Cj . 

We are treating a constant volume combustion process 
and the chemical reaction studied here ia a single 
step, irreversible reaction. 

Cj ♦ Cj ♦ Product ♦ Heat 

The reaction rate is taken to be a function only of 
the reactant concentrations and is not temperature 
dependent. Although we lose some important physical 
features of real reacting flows with this simplifi- 
cation, we are able to study details of the affects 
of heat release on the dynamics of reacting flows. 

All simulations were run at a Reynolds number 
of 250, Peclet number of 100, and a Damkbhler number 
of 1.5. Using the low Hath number approximation 
(IMA) , three different cases have been run for dif- 
ferent values of the heat release parameter, Ca, to 
study the effects of heat release on the flov. In 
addition, computer simulations involving the solu- 
tion to the exact, fully compressible (FC) aquations 
have bten performed to examine the validity of the 
low Mach number approximation and the range of hast 


release over which the approximate equationa will 
produce accurate results. 

The validity of the low Mach number approxi- 
mation in the parameter range of the simulations 
discussed here has been tested by aolving the exact, 
fully compressible equations (equations 3a - 3f) 
for a reacting mixing layer with a Mach number of 
0.2. Only a single rollup of the fundamental mode 
was computed because of the excessive computational 
costa associated with the solution of the fully com- 
pressible equations. Vorticity contours obtained 
from aolutiona of the exact and approximate equa- 
tions are compared in Fig. 1. Although acoustic 
waves propagate throughout the flow field, they 
apparently have no influence on the vorticity 
dynamics. The pressure field, plotted in Fig. 2 
shows the evidence of acoustic waves. Theae 
acoustic waves have a frequency much greater than 
the hydrodynamic frequency, so that aa mentioned 
earlier, the hydrodynamics cannot respond to the 
acoustic fluctuations, indicating that only pres- 
sure averaged over acoustic periods ia important. 

Figs. 3 and 4 are contour plots of the product 
concentration in the reacting mixing layer for runs 
1 and 3. The only difference in these two sim- 
ulations is the amount of heat release (Ce ” 0, 

Ce " 5), all other parameters being equal. The 
development of the layer into the now widely recog- 
nized large coherent structures and the pairing 
process of these structures ia clearly evident. 

The win difference observed in these two cases is 
the amount of product formed. The effect of heat 
release has been to significantly reduce the amount 
of product produced in the layer. 'Concentration 
profiles averaged horisontally across the layer are 
shown in Fig. 5 for runs 1, 2, and 3. Here we see 
the reduction in product with the degree of heat 
release. Averaged temperature profiles are dis- 
played in Fig. 6. The thickness of the layer, 
based on the point at which either the average 
product or temperature approach free stream values, 
is aeen to remain approximately the same, with a 
possible slight decrease indicated. Because of the 
decrease in density caused by exothermic chemical 
reactions, the heat release and resulting density 
changes thus act to decrease the rate of mass 
entrainment into the layer, sinee the overall layer 
growth is approximately unchanged from the cold 
flow cate. Thia is in qualitative aerement with 
the experimental findings of Wallace^®, who con- 
ducted an experiaiental study of mixing layers that 
includad a chemical reaction with weak hast release. 

Fig. 7 shows the computed velocity profiles 
across the mixing layer for runs 1, 2, and 3 at a 
time of t“24. The layer thickness, as defined by 
the point at which the average velocity attains one 
half of the free stream value, is seen to decrease 
with increasing heat release. The interesting 
feature apparent in theae plots ia the inflection 
point in the profile for the runs with heat release. 
This can ba explained in terms of the vorticity 
dynamics. For a two-dimensional mixing layer the 
vorticity aquation may be written as 

gj- * wVv - (Vp x V P ) ♦ Vx(V^) . 

P 

The first term on the right hand aide ia referred 
to as the baroclinic torque and the second term on 
the left hand side represents a change in vorticity 
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r*«ulting from thermal *xp«n»ion. In the sisnila- 
tiona dif cutted here the action of the baroclinic 
torque vat the aain source of vorticity generation. 
If there are no heat aourcea present in the flow, 

V*v « 0 and p ■ constant ao that aside from dif- 
fusion effects, vorticity is conserved along stream- 
lines. Fig. B shows the vorticity contours at timet 
of t - 8, 16, and 24 for a reacting mixing layer 
without heat release. The vorticity it of the same 
aign throughout the flow field and the large scale 
atruetures of the flow are clearly apparent. The 
vorticity field for a mixing layer with heat release 
(run 3, Ce * 5) is shown in Fig. 9. The heat 
release has resulted in the generation of vorticity 
which changes aign at the outer edges of the layer. 
This indicates a change in the sign of the velocity 
gradient, i.e. , an inflection point. The generation 
of vorticity is caused mainly by the action of the 
baroclinic torque term in the vorticity equation. 
This we have plotted in Fig. 10. A qualitative 
explanation for this complicated looking plot is as 
follows. The direction of the pressure gradient is 
roughly radially outward from the center of the vor- 
tex structures, and is monotonic across the reaction 
front. The density gradient, however, changes sign 
across the reaction front, resulting in a change in 
aign of the baroclinic torque in the reaction tone. 
As the layer rolls up, the reaction tone becomes 
highly distorted, and the complicated structure of 
the baroclinic torque shown in Fig. 10 resultst. 

Another effect of heat release in reacting flows 
will generally be to raise the viscosity, thus 
lowering the Reynolds number. It can be argued that 
this will tend to have a stabilising effect on the 
flow. If the Reynolds number is large enough the 
dynamics of the mixing layer, however, do not depend 
strongly on the Reynolds number. In the simulations 
presented here the viscosity hat been held constant. 
In other flow configurations, where the instability 
or other important features of the flow are more 
dependent on Reynolds number (for example, boundary 
layer flows), ignoring Reynolds number dependence 
could give quite misleading results. 


S. Conclusions 

The application of the approximate set of aqua- 
tions derived here for low Mach number flows has 
resulted in a significant reduction in computer 
costs compared with solutions to the exact, fully 
compressible equations of motion. The computer time 
needed to integrate the governing equations has been 
reduced by a factor of approximately 1/M using the 
low Mach number approximation. Excellent agreement 
with the exact equations has been obtained for the 
cases prasented in thia paper. 

Preliminary simulation results in two spatial 
dimensions of a chemically reacting mixing layer 
indicate thst the amount of mass entrained into the 
layer is reduced as the rate of heat release in- 
creases. The thickness of the layer is also reduced 
when exothermic chemical reactions occur. This is 
in qualitative agreement with mixing layer experi- 
ments with low rates of heat release. 


+ This qualitative picture benefited from a dis- 
cussion with Professor F. A. Williams of Princeton 
University. 


These simulations demonstrate the effectiveness 
of the direct numerical simulation approach as a 
predictive tool in studying complex reacting flows 
with heat release. An extension of this work to 
three spatial dimensions should provide insight 
into the more complex problem of the interactions 
between turbulence and heat release in reacting 
flows. 
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Abstract 


The mechanisms by which heat release affects the 
fluid dynamics in a turbulent reacting mixing layer 
are studied by direct numerical simulation. In 
agreement with previous laboratory experiments, the 
heat release is observed to lower the rate at which 
the mixing layer grows and to reduce the rate at 
which chemical products are formed. The baroclinic 
torque and thermal expansion in the mixing layer 
are shown to produce changes in the flame vortex 
structure that act to produce more diffuse vortices 
than in the constant density case, resulting in 
lower rotation rates of fluid elements. Previously 
unexplained anomalies observed in the mean velocity 
profiles of reacting jets and mixing layers are 
shown to result from vorticity generation by 
baroclinic torques. The density reductions also 
lower the generation rates of turbulent kinetic 
energy and the turbulent shear stresses, resulting 
in less turbulent mixing of fluid elements. 


changes. If the Damkohler number of the flow is 
large (fast chemistry), the reaction rate will be 
controlled by molecular diffusion and fluid 
dynamical mixing. Properly treating the turbulent 
behavior is then clearly an essential part of any 
predictive method for this type of flow. In a 
mixing layer, for example, product formation is 
augmented by stretching and wrinkling of the reac- 
tion rone. A highly strained flow field develops, 
which increases the area of the reaction interface 
and produces an inflow of reactants to this inter- 
face. For reactions that are accompanied by large 
amounts of heat release, the fluid dynamics of the 
flow will be coupled to the chemistry through the 
nonhomogeneous density field that results. The 
main objective of this work is to investigate the 
effects of exothermic chemical reactions on the 
turbulent flow field and to explain these observa- 
tions by studying the basic physical mechanisms 
that are involved. 


Calculations of the energy in the various wave 
numbers shows that the heat release has a stabil- 
izing effect on the growth rates of individual 
modes. A linear stability analysis of a simplified 
model problem confirms this, showing that low 
density fluid in the mixing region will result in a 
shift of the frequency of the unstable modes to 
lower wave numbers (longer wavelengths). The 
growth rates of the unstable modes decrease, 
contributing to the slower growth of the mixing 
layer. 


1 . Background 

The interactions between chemical heat release and 
fluid dynamics are one of the least understood 
aspects of chemicallyreacting flows. To increase 
our predictive capability of reacting flows and 
develop more reliable and efficient combustion 
models there is a need to understand more fully the 
fundamental physical processes occurring in such 
flows. In the work presented here these processes 
are studied in a temporally growing mixing layer by 
the technique of direct numerical simulation. 

Many chemically reacting flows are turbulent in 
nature and are characterized by large amounts of 
energy relesse, resulting in Isrge density 
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To achieve this purpose., direct numerical simula- 
tions of two- and three-dimensional chemically 
reacting mixing layers have been performed. An 
exothermic chemical reaction between initially 
unmixed chemical species is included. Numerical 
simulations can supplement laboratory experiments 
well and have the advantage that they can give much 
more information about the flew, since the entire 
flow field is known at every time step. In 
particular, properties that are difficult or 
impossible to measure experimentally can easily be 
obtained and studied in the simulations. In 
addition, initial conditions are easily controlled, 
and flow field parameters can easily be varied to 
study their effect on the flow. Unfortunately, 
computer time and storage requirements limit full 
simulation to flows with moderate Reynolds and 
Damkohler numbers. 

A mixing layer is formed when two initially 
separated parallel flowing fluids of different 
velocities come into contact and mix (figure 1). 
This fairly simple free shear flow has been 
extensively studied and has proven very useful in 
understanding the nature of turbulent flows. Lab- 
oratory experiments have shown that, at least in 
its early stages, the mixing lsyer is dominated by 
large-scale, quasi-two-dimensional vortices, 3 
with the growth of the mixing layer dominated by 
the pairing of these vortices . * TVro-dimensional 
simulations have been shown to accurately portray 
the characteristic large scale-rollup and vortex- 
pairing process. 3 Implications of heat releasing 
chemistry on this process can thus be meaningfully 
addressed with two-dimensional simulations. 
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However, since all turbulent flows are inherently 
three-dimensional, some potentially important 
physics will be lost by restricting the simulations 
to two spatial dimensions. For example, secondary 
instabilities can develop ir>nn streanwise vortices 
or "ribs",^ 1 which can inc.ea. fixing and enhance 
product formation. Three-dimensional simulations 
are necessary to 6tudy this type of behavior and to 
realistically treat the turbulent behavior of the 
flow. 


A large amount of work in the area of turbulent 
reacting flows has been directed at trying to 
understand the effects of the fluid dynamics on 
mixing and on reaction rates. A number of 
experiments have been performed on chemically- 
reacting constant density mixing layers to study 
the process of mixing and entrainment.®”® In 
these experiments the chemistry had no influence on 
the development of the velocity field due to the 
small amount of heat released. The results of 
these studies consistently showed that chemical 
reaction products were concentrated in large, 
spanwise-coherent structures that, at least ini- 
tially, dominate the turbulent transport. In 
addition, three-dimensional effects are present and 
can be important. Breidenthal® notes a signifi- 
cant increase in product formation rate coinciding 
with the development of three-dimensional motions 
in the flow. 


obta ined . 

Higher heat release experiments have been performed 
by Keller and Daily. ^ Cold premixed reactants 
carried in a high velocity stream were ignited by 
hot, low density combustion products which vere 
carried in a low speed stream. Results were 
reported for a range of equivalence ratios. In 
these experiments a large, favorable pressure 
gradient existed in the atreawise direction and 
the reaction was not diffusion limited. The mixing 
layer thickness was observed to increase with 
increasing heat release, a result different from 
what Hermanson and Wallace reported in their 
experiments with a diffusion limited reaction and 
an initially uniform freestream density. The 
thickening of the layer with heat release was 
attributed to downstream acceleration of the low 
speed, low density fluid. Large scale, two-dimen- 
sional vortices were again observed to persist over 
all heat release ranges. 

Two-dimensional numerical simulations of combustion 
at a rearward facing step performed by Hsiao et 
al.*® and Ghoniem et al.*® gave qualitative 
agreement with similar experiments* 7 . 

Chorin's*® random vortex method was used and 
expansion effects due to heat release were modeled 
by a set of distributed volume sources. 


Previous direct numerical simulations of an 
incompressible reacting mixing layer*® focussed 
on how the turbulent flow field affects the 
chemical reaction rate. This work has given 
physical insight into how vortex rollup and 
three-dimensional turbulence affect the chemical 
reaction. Furthermore, comparisons of these 
simulations with laboratory data have led to an 
increased confidence in the application of the 
direct numerical simulation methodology to this 
class of problems . 

Effects of heat release in turbulent reacting 
mixing layers have been studied experimentally by 
Wallace** and Hermanson.* 7 Wallace utilized a 
nitric oxide - ozone reaction to attain temperature 
rises of 400°K. Hermanson used a hydrogen - fluorine 
reaction to produce temperature rises to 940°K. 
Large-scale structures were observed to persist 
over these temperature ranges. The results that 
were obtained in these experiments all indicated 
that the heat release resulted in a slower growth 
rate of the layer and a decrease in the total 
amount of mass entrained into the layer. Hermanson 
suggested that the decrease in mass entrainment 
could be attributed to a reduction in the turbulent 
shear stresses that was observed in the high heat 
release experiments. In both sets of experiments 
the streetwise pressure gradients were small. 
Hermanson performed additional experiments while 
imposing a favorable streetwise pressure gradient. 
These experiments showed an additional thinning of 
the layer. 

McMurtry et al.*® performed direct numerical 
simulations of a two-dimensional temporally growing 
mixing layer that included effects of chemical heat 
release. A low Mach nuttier approximation was used 
which allowed them to study the effects of density 
changes vdiile eliminating the numerical stability 
requirements for computing the high frequency 
acoustic waves. Qualitative agreement with the 
experimental results of Wallace and Hermanson was 


Visual studies of diffusion flames in jets*® -7 * 
have shown that the existence of flames (heat 
release) retards the transition from laminar to 
turbulent flows. More recent turbulence measure- 
ments 77 * 7 ® show lower turbulence levels near the 
jet exit in jets with flames and an increase in 
turbulence downstream. Velocity profiles are 
steeper in the heat release runs and Yule 7 ® 
reports humps (i.e., more than one inflection 
point) in the profiles that are not seen in the 
cold jets. Also, the frequencies of the most 
energetic instabilities (vortices at the fuel-air 
interface) were observed to decrease as heat 
release increased. 

Some theoretical work performed by Marble and 
Broadwell 7 * and Karagozian 7 ® is helpful in inter- 
preting the results of our simulations. Marble and 
Broadwell studied the deformation of a constant 
density diffusion flame by turbulent motions. 
Karagozian examined the deformation of laminar 
diffusion flames in the flow field of two- and 
three-dimensional vortex structures, and also 
studied the effects of heat release on a laminar 
diffusion flame interacting with an inviscid 
vortex. The presence of density changes in the 
core caused the entire flow field to be shifted 
radially outward. The braids, or "flame arms" of 
the vortex were thus translated to a region of 
lower total straining (further from the point where 
the vortex was imposed), reducing the rate at which 
reactants were consumed by the flame. The effect 
of the reduction of the density in the core of the 
vortex structures was shown to reduce the rate at 
which reactants are consumed by the flame. 

These previous investigations have shown that there 
can be a very significant influence of combustion 
on the velocity field in reacting flows. The 
objective of the work presented in this paper is to 
use the simulation results to understand the mech- 
anisms that produce these observed heat release 
effects . 
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The Temporally Developing MixinR Layer 

The simulation results discussed here are of a 
temporally growing mixing layer which is taken to 
be statistically homogeneous in the strea. '.-f. (x) 
and spanwise (z) directions. This is not the same 
flow as the spatially developing layer that is 
usually studied experimentally, but is an approx- 
imation if one follows the flow at the mean 
velocity. By studying a temporally growing layer, 
the requirement of specifying inflow-outflow 
boundary conditions, which are difficult to 
correctly implement for the spatial layer, is 
avoided. In this work periodic boundary conditions 
are used in the homogeneous directions and pseudo- 
spectral numerical methods are used to solve the 
governing equations of motion. (See Gottlieb and 
Orszag® 2 for a description of these methods.) 

There are many dynamical features conmon to the two 
mixing layers, 26 so that a study of a temporal 
mixing layer can reveal important information about 
the spatial layer. Some of the differences between 
the two have been pointed out by Corcos and 
Sherman. In the spatially developing layer, 
events that occur downstream can induce changes in 
the flowfield upstream, whereas in the temporally 
developing layer, obviously no event can affect the 
flow at previous times. Also, the spatially 
developing layer has no symmetries around any 
spanwise axis so that entrainment rates of fluid 
into the layer from the two streams will not 
necessarily be the sane. These points should be 
kept in mind when comparing experiments to simula- 
tions. For co-flowing mixing layers the differ- 
ences between the two become small ldien the 
parameter (-(U1-U2) /(U1+U2) becomes small. U1 
and U2 are the velocities on the high- and low- 
speed sides of the layer (figure 1). To address 
these differences, simulations of a spatially 
growing layer should be performed. 

In the following section, a summary of the numer- 
ical approach used is discussed. Selected results 
that illustrate some of the effects of heat release 
on the flow are presented in section 3. These 
effects are explained in section A by studying 
different aspects of the flow including the turbu- 
lent shear stress and turbulent kinetic energy 
distribution, the vorticity dynamics, and stability 
considerations. Section 5 aummarizes the results 
and includes a discussion of how the different 
aspects of the flow discussed in section A are 
related . 


2. Methodology 

Performing numerical studies specifically directed 
at investigating the basic interactions between the 
chemistry and the flow field calls for an approach 
in tftiich the basic physical processes are an 
inherent part of the methodology. The method used 
in this work is termed direct numerical simulation 
and involves solving the three-dimensional, time- 
dependent governing equations for the detailed 
development of the flow field. This method of 
studying turbulent flows was developed by Orszag and 
Patterson 2 ® and first applied to homogeneous, 
isotropic turbulence. This technique uses no 
closure modeling and no assumptions pertaining to 
the turbulent behavior of the fluid are made. 
However, due to finite computer resources, the 
range of temporal and spatial scales that can 
accurately be resolved is limited. This limitation 


restricts these simulations to flows with moderate 
Reynolds and Damkohler nun&ers, although a con- 
served scalar approach can be used 2 ® to treat the 
limit of infinite Damkohler number. At the present 
time, direct numerical simulations have been 
limited to fairly simple geometries and are used 
only in research applications. These applica- 
tions mainly consist of trying to understand the 
physics of simple flows with the idea that the 
information obtained can then be applied to 
predicting the behavior of more complicated flows. 

Direct numerical simulations have been used 
successfully in many fluid mechanical applications 
including homogeneous turbulence, 2 ®”® 2 turbulent 
channel flow,® 2 * mixing layers.® reacting mixing 
layers without heat release, 2 ® turbulent wakes,® 2 *®® 
and turbulent boundary layers.®® Extending the 
use of direct numerical simulations to reacting 
flows with chemical heat release has been justified 
by two-dimensional simulations performed by 
McMurtry et al. 2 ® In this paper, results of 
three-dimensional simulations are presented to 
study the mechanisms by which heat release affects 
the flow. 


Low Mach Humber Approximation 


To solve for the development of the flow field, a 
set of approximate equations valid for low Mach 
number flows is utilized. This approximation has 
previously been presented, in various forms, by 
Rehm and Baum,®® Sivashinsky ,® 2 Buckmaster ,®® 

Majda and Sethian,®® and McMurtry et al. 2 ® The 
general idea behind this approximation i6 that for 
low Mach nucfcer flows, the acoustic waves generated 
by the turbulence and the combustion process have a 
much higher frequency and much faster propagation 
velocity than the motions characterizing convection 
processes. In the asymptotic limit, the speed of 
sound becomes infinite and any disturbances in 
thermodynamic pressure will be felt instantaneously 
throughout the fluid. This approximation allows 
one to study effects of density changes due to heat 
release while avoiding the numerical stability 
problems of resolving acoustic wave propagation. 

Starting with the exact equations of combustion gas 
dynamics and applying a low Mach number expansion, 
a set of ordered approximate equations is obtained. 
In this work, the equations solved include conser- 
vation equations for mats, momentum, energy, 
chemical species, and an equation of state. In the 
limit of low Mach nunfcer flows, these equations 
take the following (nondimens ional ) form: 
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The momentum equation, to the lowest order (eq. lb) 
limply describes the variation of the thermodynamic 
pressure p ^ • To complete the description of 
the velocity field the first order momentum equa- 
tion must also be retained, and is given by 


(o) d! 0) -(°). 

w Dt 
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to minimize these errors by using small time steps 
and transport coefficients (viscosity, molecular 
diffusion coefficient) that are high enough to 
ensure accurate resolution. The lack of phase and 
diffusion errors is a very desirable property in 
simulations of reacting flows, tdiere steep grad- 
ients of reacting species can develop and where 
reaction rates are controlled by molecular 
diffusion. However, truncation errors can become 
significant if gradients become too steep. 


The nondimens ional parameters appearing in these 
equations are the Danfcohler nunber, Da “ k 0 CJ. 0 /V 0 , 
Peclet number, Pe * LJJ 0 /D 0 , and Reynolds nuaber. 

Re * UoL c /v. The parameter Ce , a nondimens ional 
number characterizing the amount of heat release, 
is given by Ce ■ C a BM p 0 Cv’S 0 . 

Dq> V, C 0 , p, and T 0 are the free stream molecular 
diffusivity, kinematic viscosity, chemical species 
concentration, density, and temperature. U Q is 
the velocity difference across the layer, AH is 
the heat of reaction, and Cy is the specific heat 
at constant volume. For computational convenience, 
the length scale Iq is chosen such that the non- 
dimensional wavelength of the most unstable mode is 
2n ( L 0 ■ V2lt, where X is the dimensional 
wavelength). 

The distinction between the two pressures that 
appear, p(^) and p^ ^ , ia essential both from a 
theoretical point of view and in the numerical 
solution procedure. The thermodynamic pressure, 
p(0), is constant in space, but may vary in time 
due to heat addition, pd ) is the dynamic pressure 
associated with the fluid motion and does not 
participate directly in thermodynamic processes. 

These equations have previously been successfully 
applied to a two-dimensional problem. ^ Compar- 
isons of simulations using the exact equations and 
the low Mach number approximation were performed 
for a flow with a Mach number of 0.2 and confirmed 
the validity of the approximation in studying this 
type of flow. Details of the derivation of the 
equations and the numerical solution procedure are 
also discussed in reference 13. 

Boundary Conditions and Numerical Methods 

Since a temporally developing mixing layer is 
studied here, periodic boundary conditions can be 
applied in the streanuise (x) and spanwise (z) 
direction. In the transverse (y) direction, a 
free-slip, adiabatic boundary condition is 
applied. With these boundary conditions, vary 
accurate pseudo-spectral numerical methods can be 
efficiently implemented. All dependent variables 
are expanded in Fourier Series in the streetwise 
and spanwise directions. In the transverse 
direction a Fourier cosine series is used for all 
dependent variables except for the transverse 
velocity component, v, which is expanded in a 
Fourier sine series. Spatial derivatives are then 
computed by differentiating the series term by 
term. A second order accurate Adams-Bashforth 
time-stepping scheme it used to advance the 
equations in time. 


Flow Field Initialization 

The computational domain for the simulations is 
chosen to be large enough to contain the most 
unstable mode and its subharmonic (as determined 
from linear stability theory for an incompressible 
flow^O). The velocity field is initialized by 
adding a hyperbolic tangent velocity profile 
(figure 2) to a low level, three-dimensional, broad 
band background perturbation. In addition, in 
certain simulations a perturbation in the form of 
the most unstable mode and its subharmonic are 
added. The spatially periodic "forcing" in the 
temporally developing layer is analogous to period- 
ic (in time) forcing of a spatially growing layer. 

To specify the background perturbation we use a 
method similar to that discussed by Orszag and 
Pao. 22 (Also discussed in more detail in Riley 
and Metcalfe.33) Energy is assigned to each wave 
number component as specified by a selected three- 
dimensional energy spectrum. A random phase is then 
given to each component which results in a random 
velocity field with a specified energy spectrum. 

The actual shape of the spectrum is not critical as 
long as energy is contained in a fairly wide range 
of wave numbers. The spectrum used in these 
simulations is given by 

4 ,4 

E^k) “ CjjA 2 2 3 
( l+k^A 2 ) 

This is a fairly broad banded, isotropic spectrum. 

A is the integral length scale, k is the magnitude 
of the wave number, and £33 is a coefficient 
that determines the level of the perturbation. 

This spectrum is multiplied by a "form function" in 
physical space to give a turbulence intensity 
profile characteristic of mixing layers. ^ For 
simulations that start with a low enough initial 
amplitude (in the linear range), the most amplified 
frequencies will ultimately dominate. 

The chemical reactant fields are initially 
two-dimensional and are given by the following 
functional relationships: 

y 

C^.O) /expf-C^dC 

•<8 

C 2 (x,y,*,0) - 1 - C 1 (x,y+6,z,0) 


The pseudo-spectral technique as implemented here 
exhibits very small phase errors and numerical 
diffusion compared to finite difference techniques. 
The only errors that are of any consequence are 
truncation errors due to the finite wave nunber cut 
off and time stepping errors. Care has been taken 


£ is a value by which the two species fields have 
been offset so that there is initially no overlap 
betveen the two (figure 2). These functions are 
easily resolvable using spectral methods and are 
not unlike those measured experimentally. The 
particular chemical reaction used is a single step. 


39 



irreversible reaction. 


Cl * Cj ■* Products 4 Heat 

The reaction rate is a function only of the 
reactant concentrations and does not depend on 
temperature. Although some important features of 
real reacting flows are lost with thi6 simple 
reaction, many important features of the effects of 
energy release on the dynamics of the flow can 
still be studied. 

Summary of Approximations 

The following is a summary of the major simpli- 
fications and approximations used in the numerical 
simulations . 

A. low Mach number approximation. 

To ease stability requirements in the numerical 
calculation a low Mach number approximation is 
used. This allows acoustic waves to be filtered 
out of the equations. The validity of the 
approximation has been illustrated for these 
simulations in reference 13 by comparing results 
from the complete set of full compressible 
equations with the results from the low Mach number 
approximation . 

B. Temporally developing mixing layer. 

The simulations are all for a mixing layer that 
develops in time and is spatially periodic. As 
pointed out earlier, this is not the same flow as 
the spatially growing layer that is encountered in 
practice. However, many dynamic similarities exist 
between the two so the temporally growing case is 
of relevence. Also, a more efficient and accurate 
computer code can be written for the temporally 
growing mixing layer for a given amount of computer 
res ources . 

C. Simple chemistry. 

A single step, irreversible reaction between two 
species to form a product is used. The reaction it 
taken to be a function only of the reactant 
concentrations and does not depend on temperature . 
This does not represent any real chemical reaction 
but does allow for details of the effects of energy 
release in the flow to be studied. 

D. Constant transport coefficients. 

The viscosity, thermal and molecular diffusivities , 
and the specific heats are taken to be temperature 
independent constants. Again, this is not a 
realistic feature of general cocbusting flows, but 
in addition to simplifying matters from a numerical 
point of view, it allows other effects on the flow 
to be isolated and studied in a simpler environment.- 

E. The amount of heat release is restricted. 

This restriction is necessary to ensure that large 
velocities will not be generated by the cozfcustion 
that may violate the lew Mach number 
approximation. The condition that must be 
statisfied is DaCe 1/M . 


3. Simulation Results 

Results are reported for two different values of 
the heat release parameter, Ce: one with no heat 
release (Ce“0), and the other giving a density 
decresse of p/p 0 - 0.5. Runs 1 (Ce«0) and 2 
(Ce"5) were initialized with random perturbations 


only. Runs 3 (Ce"0) ana 4 (Ce«5) included, in 
addition, a two-dimensional perturbation 
corresponding to the most unstable wavelength and 
its subharmonic. In the follcving, runs 3 and 4 
will be refer-eu <o ls the "forced" cases. These 
runs are analogous to laboratory experiments in 
which well defined harmonic perturbations are 
applied to the flow at selected frequencies. The 
parameters used in these runs are given in table 
1. Additional simulations have been run using 
different random initializations. The results 
presented here are representative of all these runs. 

All simulations were carried out on the CRAY X-MP 
computer at NASA Lewis Research Center. The three- 
dimensional computations were performed on a 64 x 
65 x 64 array and required 12 seconds of epu time 
per time step. Each simulation presented here 
required between 600 and 1200 time steps to 
complete . 

The initial three-dimensional random perturubations 
added to the velocity field were identical for all 
three-dimensional runs. The initial energy spec- 
trum for runs 1 and 2 (three-dimensional random 
perturbations only) and runs 3 and 4 ( which 
includes the two-dimensional forcing) is shown in 
figure 3. The energy contained in the fundamental 
wavelength (wavenumber ■ 1 ) is approximately an 
order of magnitude higher in the forced cases, 
although the energy in the higher modes is the 
same. The initial streamwise rms velocity profile 
is shown in figure 4 for the forced runs (run 3 and 
4). The initial turbulence levels for these runs 
is low, approximately 32 of the velocity difference 
across the layer. Tor the unforced runs (runs 1 and 
2), the maximum turbulence levels were 
approximately 2.22 . 

Higher heat release cases can be run, but since the 
boundary conditions used here imply constant volume 
combustion, the equivalent of a significant down- 
stream pressure gradient would develop (the back- 
ground pressure would rise with time). Test cases 
with a computational domain twice as large in the 
transverse (y) direction (reducing the effective 
pressure gradient by a factor of two) were perform- 
ed to verify that pressure gradient effects were 
not significant for the values of the paramters 
s tudied here. 

The grevth of mixing layers is strongly influenced 
by large-scale two-dimensional vortices , and to 
understand the development of the flow field, the 
dynamics of these structures must be understood. 
Figure 5 is a contour plot of the instantaneous 
• panwise vorticity field obtained from a previous, 
constant density, two-dimensional simulation*^. 

The development of the most unstable mode and the 
pairing of two vortices to form a single larger one 
is apparent. In this simulation the velocity field 
was initialized with perturbations corresponding to 
the most unstable mode and its subharmonic. One 
purpose of this present investigation is to examine 
the influence of heat release on this process and 
the influence this has on the product formation 
rate. 

Some aspects of the three-dimensional nature of the 
flow can be seen in a three-dimensional perspective 
plot of the total vorticity. In figure 6, surfaces 
having a value of 502 of the maximum of the sum of 
the absolute values of all three vorticity compo- 
nents are plotted ( lo^l 4 luyl 4 I ■ constant). 
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This is a result from run 1 (initial random velo- 
city perturbations, no heat release) after one 
pairing. In addition to the approximately two- 
dimensional apanwise vorticity associated with the 
large-scale structure, tubes aligned in the 
streaswise direction are also apparent. These 
structures have been recognized as counter-rotating 
vortices 4 and have been modeled by Lin and 
Corcos^® and computed in constant density mixing 
layers by Metcalfe et si. 26 The effect of these 
structures on the flow field it to induce a 
velocity field that acts to pump fluid between the 
two layers, thereby further convoluting the 
reaction interface. This tends to increase mixing 
between the two streams and enhance chemical 
reactions. 

The effect of heat release on the chemical product 
formation for runs 1 and 2 ia shown in figure 7, 
where the total product (integrated over the 
entire computational domain) is plotted as a 
function of time. In agreement with previoua 
two-dimensional results, 1 ^ the total product ia 
seen to be lower for the case with heat release. 
Since the instantaneous reaction rate (the slope of 
these curves) is only a function of the joint 
species concentrations, this again indicates that 
one effect of the heat release has been to reduce 
the amount of mixing between the two streams. An 
increase in the rate of product formation is seen 
to occur at about t • 40 in the constant density 
case (run 1) and at t ■ 60 in the heat release 
case. This is due to the rollup of the subharmonic 
mode, which engulfs large amounts of fluid into the 
vortex structures. Figure 8 shows the amount of 
product formed when the layer is forced at its most 
unstable wavelength and subharmonic in addition to 
the random three-dimensional perturbations. The 
overall rate of product formation is greater when 
the layer is forced, but again the decrease in the 
product with heat release is apparent. 

The total product from a two-dimensional ainulation 
with the same level of forcing as runs 3 and 4 
(Aj o * 0.01 and * 1 / 2,0 “ 0.01) is also plotted 
in figure 8. During the rollup and pairing pro- 
cess, the three-dimensional growth ia inhibited, 
giving nearly identical results for the two- and 
three-dimensional simulations. Other studies 2 ^ 
have shewn rollup and pairing has a stabilizing 
effect on the three-dimensional growth when the 
amplitudes of the three-dimensional modes are 
small, vdiile absence of pairing can enhance the 
grewth of the three-dimensional modes. After 
saturation of the two-dimensional modes (at a time 
of about t«35) the product formation in the 
three-dimensional runs continues to grow rapidly 
due to enhanced three-dimensional mixing, while in 
the two-dimensional simulations the product 
formation drops off rapidly. 

One measure of the layer growth ia the vorticity 
thickness, defined as the inverse of the maxisum 
slope of the velocity profile. The velocity 
profiles for runt 1 and 2 (figure 9) at a time of 
t<*72 and runs 3 and 4 (figure 10) at a time of t"36 
show a steeper profile for the heat release runs, 
indicating a slower rate of the layer grwth. Note 
also the overshoot in the velocity profile that 
occurs in the heat release runs. These observa- 
tions are similar to those obtained from previous 
two-dimensional calculations, 11 except the 
overshoot in the velocity profile is not as 
pronounced in the unforced three-dimensional 


simulations (runs 1 and 2) as in the simulations 
with two-dimensional forcing. This is because of 
the greater apanwise variation and lack of 
coherence that exists in the three-dimensional 
s inulations . These overshoots seen in the velocity 
profile are similar to those sac-, iii the jet flame 
experiments of Yule et al. 2j and have also been 
Been in reacting mixing layer experiments. 48 

Plots of the vorticity thickness and velocity 
half-width (the location where the average stream- 
wise velocity component attains one half of its 
free stream value) as a function of time reveal 
another interesting feature (figures 11 - 14). 
Initially, the grewth of the layer given by these 
two measures is slightly greater in the heat 
release case; up to t a 30 in the unforced runs (runs 
1 and 2, figures 11 and 12) and up to t“15 in the 
forced runs (runs 3 and 4, figures 13 and 14). At 
later times the layer thickness is less in the heat 
release runs. In addition, the initial difference 
in the grewth rates with heat release is not as 
great in the forced runs as in the unforced runs. 
The initial increase in the layer growth rate is a 
result of thermal expansion, sAiich shifts the uhole 
flowfield outward. Explanations for the smaller 
thickness seen in the heat release runs at the 
later times, as well as the differences seen in the 
forced and unforced cases are given in the next 
section. 

To summarize, the roost obvious macroscopic effects 
of heat release on the mixing layer dynamics 
revealed in these numerical simulations and similar 
experiments are a decrease in the product formation 
and, after a small initial increase, a decrease in 
the layer growth rate. These results are 
qualitatively similar to those obtained experiment- 
ally by Wallace 11 and Hermanson, 12 and are 
consequences of lower rates of fluid entrainment 
into the mixing region when exothermic chemical 
reactions occur. In the following section, physical 
mechanisms responsible for these observations are 
made . 


4. Effects of Heat Release 

In this section, the evolution of the flow and the 
influence of heat release is discussed in terms of 
the shear stress distribution, the turbulent 
kinetic energy balance, vorticity dynamics, and 
stability considerations. Since each of these 
aspects of the flow must be consistent with the 
others, they are not interpreted as separate 
mechanisms, but instead they provide different 
viewpoints from which the effects of heat release 
can be explained. 

Turbulent Stresses 

Hermanson 12 suggests that most of the observed 
heat release effects can be accounted for as a 
result of a decrease in the turbulent sheer 
stresses. Velocity and density profiles obtained 
experimentally as well as the measured layer grcxvth 
rate were used to compute the turbulent shear 
stress profiles in his experiments. 

To understand the significance of the turbulent 
stress distribution it is instructive to look at 
the averaged momentum conservation equation: 
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The above equation is written in its Favre averaged 
(density weighted) form. The wavy overbar denotes 
Favre averaged quantities, which are defined as 


U - pu/p 

The fluctuating quantities are then given by 
u" " u - U . 


The turbulent kinetic energy profile (pi'.V.') for 
runs with and without heat release it shown in figure 
17 at a nonditnensional time of t“48. (Time jr. 
nondimens ionalized by L 0 /U 0 .) This corresponds 
to a time after the rollup of the most unstable 
mode and before the rollup of its subharmonic is 
complete. From this figure it can be seen that the 
total turbulent kinetic energy is less for the heat 
release run. This is consistent with the earlier 
observations of less product formation and lower 
growth rates, since lower turbulence levels imply 
lower (turbulent) transport rates, resulting in a 
lower exchange of mass, momentum, and energy among 
fluid elements. 


The straight overbar refers to conventional Reynolds 
averaging. In the temporal calculations this is sn 
average overa horizontal plane. The turbulent 
atresses, pjJVTV , which appear in the averaged 
momentum equation represent a transport of mean 
momentum (per unit volume). The Favre averaged 
atresses contain momentum exchange mechanisms due 
to turbulent transport, interactions between the 
mean flow and volume trie changes, and fluctuation- 

fluctuation interactions (p'ujut).* 1 * In the 
numerical simulations performed here , it was 
possible to compute the Favre averaged turbulent 
shear stresses directly. 

The computed Favre (density we ighteji), aver aged 
turbulent shear stress profiles (pjVulj) for runs 
with and without heat release are shown in figures 
15 and 16. Figure 15 is for a three-dimensional 
calculation without forcing (runs 1 and 2) and 
figure 16 is the profile computed for a three 
dimensional calculation that includes forcing at 
the most unstable mode and its subharmonic (runs 3 
and 4). The suppression of the shear stress in the 
heat release runs is clearly indicated. 

The lower turbulent stresses in the heat release 
cate imply a lower transport of momentum to the 
turbulence. In the following sections it will be 
shown that this also indicates that less energy is 
being transferred from the mean flow to the turbu- 
lence. In addition, the turbulent shear stresses 
can be directly related to the stability character- 
istics of the mixing layer, providing explanations 
for the lower growth rates. 

Turbulent Kinetic Energy 


Another aspect of the flow that is useful to 
examine and is also closely related to the turbu- 
lent shear stresses 'is the turbulent kinetic 
energy. This can give another way to interpret the 
effects of heat release on the turbulent motions 
and account for some of the observations in the 
heat release runs. Furthermore, the turbulent 
kinetic energy and its production, redistribution, 
and viscous dissipation are important aspects of 
the flow that must be treated in many of the models 
currently used to describe turbulent flows. 

Statistical information related to the turbulent 
kinetic energy was computed for all the simula- 
tions. Qualitatively, the heat release effects 
were the same for the cases with and without 
two-dimensional forcing. Therefore, to simplify 
the discussion, only the runs with the initial 
random velocity perturbations (runt 1 and 2) are 
discussed below. 


To understand how a lower turbulent kinetic energy 
profile results, it is useful to examine its trans- 
port equation. In the Favre averaged form, the 
transport equation for the turbulent kinetic energy 
it given by 


‘5T < * ) 4 -a^ (pu i q) 




- , ZL 
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where q * 1/2 u'.'uV is the turbulent kinetic 
energy per unit ma^s . 


( 3 ) 


From equation 3 we can see that a number of 
different mechanisms contribute to the kinetic 
energy balance. 

The first term on the right hand side is the 
production term and accounts for the exchange of 
energy between the mean flow and the fluctuating 
motion. The production of kinetic energy is 
proportional to the average velocity gradient and 
the turbulent stresses. In the mixing layer 
similated here, the only nonzero contribution of 
this term is for i«l and k”2. The profile of the 
turbulent production is shown in figure 18 at 
t“48. Production is greatest at the center of the 
mixing layer idiere the velocity gradient and 
turbulent stresses are the largest. (In the 
following plots, negative values indicate a 
production of turbulent energy). Although it was 
ahown that the velocity gradient is steeper in the 
cases with heat release, the total turbulent 
kinetic energy production is less with heat release 
due to the smaller turbulent stress term, 
pilSil’, as shewn in the previous section. The 
slope of the velocity profile is plotted in figure 
9 and the Favre averaged turbulent stress profile 
for this run is shown in figure 15. 


The next term, •— -(pu'.’u’.V') is conservative 
® x k 1 1 K ♦ 

(in general, any term of the form V*^ is conser- 
vative for a variety of boundary conditions) and 
describes the transport of turbulent kinetic energy 
by the fluctuating velocity field. From figure 19 
we see that turbulent transport tends to convect 
kinetic energy away from the center of the mixing 
layer, where the turbulent intensity is highest, to 
the outer regions. The magnitude of this term is 
less for the heat release runs since the turbulent 
fluctuations are lower. 
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The terra u'.' describes the effect* of pressure 
1 3x. 

fluctuations on the turbulent kinetic energy 
distribution. The value of this terra is plotted in 
figure 20 at t«48. This tern is of the same order 
of magnitude as the production tern. It has an 
opposite sign, however, and does not appear to be 
as greatly affected by the heat release, although a 
decrease in the magnitude is apparent. The 
physical interpretation of this contribution can be 
clarified somewhat by writing this term as 





(4) 


Another approach to studying the flow and inter- 
preting the effects of heat release on the flow 
field is in terms of vorticity dynamics. In a 
two-diraens ional flow without heat release or 
diffusion, the vorticity is conserved following 
particle paths. Therefore, observing the vorticity 
field allows a direct v isualixation of the flow 
field. In a three-dimensional flow, or a flow with 
density variations or expansion, the vorticity no 
longer serves as a reliable fluid marker. However, 
the dynamics of the flow field can still be 
understood and interpreted by studying the 
vorticity field. In this section, the effects of 
heat release on the flow field are discussed and 
explained in terms of the vorticity dynamics. 


The first terra on the right hand side is a conser- 
vative term and represents the redistribution of 
kinetic energy by pressure fluctuations . The 
second term, which is aero for constant density 
flows, is a source of kinetic energy resulting 
completely from the combustion. The contribution 
of these two terms for the heat release run is 
shorn in figure 21. The transport due to the 
pressure fluctuations is not such changed between 
the runs with and without heat release, although 
the peak magnitude is somewhat lets in the heat 
release run. In both cases pressure fluctuations 
act to transport energy from the center of the 
mixing layer to the outer regions. The contribu- 
tion due to the velocity divergence leads to a 
production of kinetic energy in the center of the 
vortices. In this respect the heat release acts to 
increase the turbulence level, although this effect 
is dominated by the decrease in the mean flow 
production term. In the outer regions, the 
pressure-divergence correlation changes sign and 
acts to decrease the turbulence level. 

From the results presented here (figures 17-21) it 
is seen that the most important effects of heat 
release are significant reductions in the turbu- 
lence production and turbulent transport. In the 
case of the production term, the lower density 
resulting from the combustion is directly respon- 
sible for most of this change. Transport by the 
fluctuating motion decreases in the heat release 
runs simply because the turbulence levels are 
lower. The heat release results in a production of 
turbulent kinetic energy along the center of the 
mixing layer through the expansion part of the 
velocity-pressure gradient correlation. However, 
the production is small in these simulations 
compared with the decrease in the exchange of 
energy with the mean flow, yielding an overall 
lower turbulent kinetic energy profile. As 
heat release is increased, the term 3u^/3xf would 
also increase , resulting in more internal energy 
being converted to kinetic energy. With 
significantly higher heat release, this could 
possibly result in an overall increase in the 
kinetic energy. 

Vorticity Dynamics 

The numerical simulations can provide us vith such 
more information about the flow than ia revealed in 
the turbulent stress profiles and turbulent kinetic 
energy profiles, which only give integrated effect* 
and do not reveal much about the physical mechan- 
isms at work. It is more instructive to relate the 
heat release effects directly to the dynamics of 
the large-scale structures. 


Vorticity Equation 

The vorticity equation is derived by taking the 
curl of the momentum equation. In a general three- 
dimensional flow, the vorticity equation can be 
written in the following form: 

“ (w»V)v - w(V*v) ♦ (Vp x Vp) / p* * -g^(V 2 u) (5) 


Four different mechanisms can be identified that 
alter the development o£ the vorticity field: 
^ortjx stretching (o*V)v, thermal expansion, 

w(V*v), baroclinic torques. (Vox Vp)/p^, and 
viscous diffusion. For a low Mach number flow 
without heat release, the expansion and baroclinic 
torque terms will be zero. When density changes 
due to heat release do occur, these two mechanisms 
can be important in the vorticity dynamics. For a 
two-dimensional flow, the vortex stretching term is 
identically zero, and if no heat ia released 
equation 5 becomes 





( 6 ) 


showing that vorticity follows fluid particle paths 
in the absence of diffusion. 


The instantaneous spanwise component of vorticity 
(U3) at times of t“48 and 72 is shown in figures 22 
and 23 for runs 1 and 2 (three-dimensional pertur- 
bations, no forcing). In the following figures, 
dashed contour lines indicate negative vorticity 
(local rotation of fluid elements in the clockwise 
direction) and aolid lines indicate positive 
vorticity. One of the most apparent differences 
between these two figures is that, for the simula- 
tions with heat release, the maximum amplitude of 
the vorticity, which occur* in the vortex cores, 
has decreased substantially. Furthermore, the 
vorticity ia not as concentrated in the center of 
the large structures as in the constant density 
case. Also note the regions of positively signed 
vorticity that appear at the outer edges of the 
vortex structures for the case with heat release. 
This same behavior has been seen in previous two- 
dimensional calculations. 1 - 3 The results of the 
cates with and without forcing are again qualita- 
tively similar so only the results with the initial 
random velocity field (runs 1 and 2) are discussed 
below. In the following, the mechanisms responsible 
for these effects and their influence on' the flow 
field development are discussed. 
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Thermal Expansion and Baroclinic Torque 

In an expanding flow, V*v it positive; therefore, 
the effect of thermal expansion results in a 
decrease in the magnitude of vorticity (see 
equation 5). This can also be understood by 
angular momentum considerations , since as a fluid 
element expands due to heat release, the magnitude 
of vorticity must decrease to conserve angular 
momentum. 

The instantaneous value of the spanwise component 
of the expansion term is shown in figure 24 at two 
different times. Thermal expansion occurs as heat 
is released by the chemical reaction so that this 
term alto gives a good indication of the reaction 
tone. The rate of reaction is highest where the 
inflow of reactants is the highest. This occurs in 
the regions of highest strain, tftiich are located in 
the braids. The result of the expansion is a 
decrease in the magnitude of vorticity, tdiich in 
part explains the changes in the vorticity field 
that result when exothermic chemical reactions 
occur. (For a temperature dependent reaction, the 
high dissipation rates in the braids can cause 
local quenching of the flame so that the reaction 
rate is not necessarily highest in the regions 
where the inflow of reactants is the highest J* 1 * 

This effect is not addressed here.) 

Another mechanism that affects the vorticity 
distribution in the mixing layer is the baroclinic 
torque, (Vp x Vp ) / c? . This results from nonaligned 
pressure gradients and density gradients. Plots of 
the spanwise component of this term at specific 
spanwise locations show an alternating sign across 
the reaction surface and no contribution in the 
vortex cores (figure 25). This can be understood 
by recognizing that the density gradient changes 
sign across the reaction surface and the pressure 
gradient vector points approximately radially 
outward from the vortex cores. In the vortex cores 
the pressure and density gradients are small and 
are also approximately aligned, to there is no 
contribution of this term here. As the layer rolls 
up and winds around itself, the baroclinic torque 
takes on the complicated structure shown in figure 
25b. Comparisons with the results of two-dimen- 
sional simulations show that these effects are 
dominated by two-dimensional dynamics. 

Comparing figures of the spanwise vorticity compo- 
nent (figure 22) and the instantaneous baroclinic 
torque (figure 25) shows that the regions of 
positively signed vorticity at the outer regions of 
the vortices, and also the appearance of multiple 
extrema in the vorticity field, are a result of the 
baroclinic torque. The overshoot in the velocity 
profile seen in the previous section (figures 9 and 
10) ia a result of the generation of positive 
vorticity in this region by the baroclinic torque. 
This is also the mechanism that produces the 
previously unexplained inflection points seen in 
the velocity profile at the outer edges of jet 
diffusion flames. 23 

In an attempt to understand the relative importance 
of the expansion and baroclinic torque terms on the 
development of the flow, average values of these 
two terms are plotted as a function of the height 
across the mixing layer for a sequence of times 
(figure 26). (Negative values of the expansion term 
indicate a decrease of vorticity while negative 
values of the baroclinic torque indicate an 


increase.) The magnitude of the expansion term is 
seen to be initially stronger than the baroclinic 
torque term. At later times, the amplitude of the 
two terms are comparable. Note that the expansion 
term consistently produces a net decrease in 
vorticity. The baroclinic torque tends to decrease 
the vorticity near the upper and lower limits of 
the dynamically active region, at least during the 
initial stages when the two-dimensional modes 
dominate the flow. This is consistent with the 
generation of regions of positive vorticity at the 
edges of the cores (figure 25). Note also that the 
baroclinic vorticity generation on the centerline 
is occurring in the braids rather than the cores. 

The baroclinic torque and thermal expansion have 
shown to result in weaker and more diffuse vortex 
structures. This will result in a slower rollup of 
the layer, thus reducing the straining of the reac- 
tion interface and decreasing the mass entrainment 
into the layer. This accounts for the decrease seen 
in the overall product formation and the changes in 
the layer growth rate. In section 3 it was shown 
(figures 11-14) that with heat release, the layer 
growth rate initially increases, and then decreases 
compared with the constant density case. The ini- 
tial increase was explained to be due to thermal 
expansion, which tends to shift the whole layer 
outward. Later in the development of the mixing 
layer, however, the growth is dominated by the 
large scale rollup and pairing process. Since this 
process is inhibited by heat release, the constant 
density mixing layer will then grow faster. Note 
also that this behavior is more apparent between 
runs 1 and 2 (unforced) than between runs 3 and 4 
(forced). In runs 3 and 4, the effects of two- 
dimensional rollup are felt sooner, since the flow- 
field was initialized with coherent perturbations 
(in addition to the random background noise) cor- 
responding to the most unstable wavelengths. The 
initial rms amplitudes of the velocity perturbations 
were approximately 25Z higher in runs 3 and 4 than 
in runs 1 and 2, due to the addition of the coherent 
two-dimensional perturbations. 

The relative importance of the expansion and the 
baroclinic torque on the flow field development 
depends on the geometry of the mixing layer. 
Experiments and simulations have shown that, in a 
reacting mixing layer, eonbustion products are 
concentrated in the vortex cores. This implies 
that the density changes in the cores will be 
responsible for most of the observed heat release 
effects. Since the baroclinic torque is weakest in 
the cores, it might be expected that thermal 
expansion alone can account for most of the 
significant changes of the altered vorticity 
field. Furthermore, from figure 25 it can be seen 
that regions of strong vorticity generation due to 
the baroclinic torque are located adjacent to 
corresponding regions of vorticity destruction. 

This is due to the change in sign of the density 
gradient across the reaction surface as explained 
above. Therefore, although the magnitude of the 
baroclinic torque may be quite large, the effect on 
the flow field falls off rapidly with distance from 
these regions. This may cause local changes in the 
flow field, but the global effects on the large 
scale structures will not be as strong as those due 
to thermal expansion. However, as the vorticity 
generated by the baroclinic torque is corivected in 
the layer and engulfed into the large scale 
structures, the integrated effects over time of the 
baroclinic torque and thermal expansion become 


44 



difficult to isolate. Since the local magnitudes 
of these two terms are of the same order both 
effects should be accounted for. 

The effects of heat release on the three- 
dimensional, spanwise variation in the reactant 
concentrations presented in section III can also be 
interpreted in terms of the vorticity dynamics. 
Secondary instabilities in the flow are character- 
ized by counterrotating streanwise vortices.^ 

The flow field induced by these vortices enhances 
mixing and increases the surface area of the 
interface between the two reacting chemical 
species, as shown in figures 27 and 28. In figure 
27 the streanwise component of vorticity at a time 
of t”72 is plotted at the streacuise location x * 
it. Comparing this with a similar plot of run 2 
(figure 28) shows a less intense streacuise 
component of vorticity tdien heat release accom- 
panies the chemical reaction. The changes seen 
between these two figures are typical of any 
streacuise cut. Undoubtedly, mechanisms similar to 
those discussed for primary instabilities also work 
to reduce the growth of secondary instabilities in 
the presence of heat release. 

Stability Considerations 

A question that arises in these simulations is 'how 
are the stability characteristics of the mixing 
layer affected by heat release?' There are both 
physical and numerical reasons for addressing this 
question. If the unstable modes shift to other 
frequencies, or if the growth rates change as a 
result of density changes, this certainly can 
affect the growth of the mixing layer and the rate 
at which chemical products are formed. 

From a computational point of view, it is not 
possible to look at a continuous distribution of 
frequencies. Because of the periodic boundary 
conditions employed here, only disturbances of 
wavelengths that divide exactly into the computa- 
tional domain are allowed (X~domain size/n). 
Therefore, if the most unstable modes shift to 
different wavelengths , dynamically important 
effects may be overlooked. 

To address this question, a linear stability analy- 
sis of a simplified model problem was carried out 
in conjunction with the simulations. This analysis 
involved a shear layer with a piece-wise linear 
velocity profile with a low density in the velocity 
transition region (figure 29). The velocity and 
density in the free stream were constant, and the 
density in the transition tone was also constant 
and given by p(l-D), 0 <33 <1 . Although this does 
not exactly describe the conditions in the aisula- 
tions reported here, the basic characteristics of 
the two are similar, so that the general trends 
indicated by the analysis are expected to be true 
of the aimulations. Details of the analysis are 
given in Me Mur try 5 

The results of this analysis arc plotted in figure 
30, which shows the growth rate of any individual 
unstable mode (specified by its wave number) for a 
given value of 0. As the density decreases 
(increasing Q), the most unstable mode, represented 
as the maximum value on each curve, shifts to a 
lower wave number (longer wavelength). The grevth 
rates of the unstable modes are also seen to 
decrease as the density is lowered, a result 
consistent with the lower growth rates discussed in 


section 3. 

The energy, E(k x ,k z ) contained in various inodes 
as a function of time is shevn in figures 31 and 32 
f oi runs 1 and 2 (initial random perturbation 
velocity field). The energies are defined here as 


E(k x ,k z ) ■ J h/(k x ,k z ,y) l 2 d) 


where N 

v (kx .kj ,y ) ■ u (k x ,k z ,ky ,t )X(ky ,y ) 

k y «0 

and u (k,t ) are the Fourier components of the 
velocity and X is either sin(k y y) or cos(k y ,y), 
depending on the component of the velocity under 
consideration. In the following discussion , we 
analyze the temporal behavior of the energies in 
these Fourier modes. Although we refer to the modes 
Ei o • nd ^1/2 0 ** the fundamental and subharmonic 
respectively, it is important to note in the inter- 
pretation of these results that this correspondence 
it not exact. For example, the subharmonic mode by 
itself does develop higher wsvenunfcer components as 
it rolls up. Thus, the Fourier mode Ej q, while 
being dominated by the energy of the fundamental , 
can also contain some energy of the subharmonic 
mode . 


In figure 31, the the fundamental (Ej q) grows 
until a time o£ t ■ 40, at which point it reaches a 
quasi -equilibrium, saturated state. 3,26 
subharnonic (X \/2 0^ continues to grew until it 
reaches its saturation level at t ■ 65. The 
behavior of these two modes changes significantly 
when heat release occurs. The growth of the 
fundamental drops off at a level well below the 
incompressible saturation level and at a much 
earlier time (figure 33), although the initial 
growth before much density change occurs (up to 
about t"7) is the same in the two cases. The 
subharmonic remains unstable but grows at a lower 
rate and reaches saturation at a later time than in 
the constant density case (figure 34). The energy 
in the three-dimensional modes (the sum of E(k x ,k z ) 
for all k x and k z f 0) is also decreased due to the 
heat release as shown in figure 35. 


The results of the simplified linear analysis for 
the linear profiles compare favorably with the 
simulation results. The density changes occurring 
in run 2 correspond to a Q of 0.5. With this 
amount of heat release, the most unstable mode in 
the constant density case is predicted to be 
stable, which is consistent with the modal behavior 
in the siculations (figure 33). The actual most 
unstable mode shifts to a wove number of 0.3, which 
is fairly close to the wave number of the 
subharmonic of the most unstable mode in the 
constant density case (0.22). From figure 34, the 
growth rate of the subharmonic is shown to reduce 
from 0.16 to 0.13, while in the simulations the 
com{uted grewth rate of the subharmonic decreased 
from 0.12 to 0.09 in the heat release case. In 
figure 34 this leaver growth rate of the subharmonic 
is indicated by the smaller slope. 


The temporal behavior of the various modes when the 
layer is initially forced at the fundamental and 
subharmonic (runs 3 and 4) is shown in figures 36 
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and 37. With these levels of forcing, the grcvth 
of the fundamental it not auppretaed by the heat 
release nearly as dramatically as in the case with 
no forcing. The fundamental stabilizes at a 
slightly earlier time in the heat release case (at 
t“18 vs. t"21 for the constant density case) and 
the energy contained in this node is only a factor 
of two lower. This observation may be a result of 
the time it takes for the density decreases to 
develop in the simulations. During this time, the 
instabilities will grcv. With high enough levels 
of forcing, substantial rollup can occur before the 
stabilizing effects of heat release are strongly 
felt. 

The laboratory experiments on mixing layers, 
although clearly showing lower growth rates as heat 
release is increased, do not conclusively show 
anything regarding the total supression of some 
modes or a shift in the wavelengths of the most 
unstable modes. This particular aspect has not 
been carefully investigated in the experiments, 
although Hermanson 12 auggests there is little 
change in the spacing (wavelengths) of the vortex 
cores when heat release occurs. This has been 
observed, however, in experiments on jets (Yule et 
al., 1981) which showed the value of the most 
energetic frequencies decreased as heat release 
increased. 

There are a number of possible reasons for why this 
behavior may not be as apparent in laboratory 
experiments as in the simulations. First, the 
suppression of higher wave number components in the 
simulations could in part result from the thicker 
reaction zone that exists in the simulations com- 
pared to the experiments. This is a consequence of 
a lever Danfcohler nunber in the simulations, which 
is necessary to achieve accurate resolution of the 
reaction zone structure. Secondly, in the exper- 
iments, the turbulence levels in the boundary layer 
at the splitter plate are often much higher than in 
the sieulations presented here, so that the flow 
instabilities may not be governed by linear stabity 
theory. Finally, any laboratory experiment on 
mixing layers is extremely sensitive to any type of 
external forcing. This forcing is often applied 
purposely, or can result from any disturbances or 
resonances associated with the experimental setup. 
Such resonances exist in any experimental facility, 
and even extremely low levels of facility generated 
coherent disturbances can result in large vari- 
ations in the development of the shear layer 
The sieulations performed here (runs 3 and 4) also 
show that even low levels of coherent forcing have 
a significant impact on the growth of the unstable 
modes. This problem of how heat release affects 
the grcvth rates and wavelengths of the unstable 
modes it an area that calls for further 
experimental work. 

5. Summary 

The simulations performed as a part of this work 
and similar laboratory experiments have shown that, 
when heat release accompanies the chemical reac- 
tion, the mixing layer grows at a lower rate and 
the amount of product formed decreases. In addi- 
tion , an over shoot in the velocity profile appears. 
In the previous section, four different aspects of 
the flow were studied to explain the observed 
effects of heat release: the turbulent shear 

stresses, the turbulent kinetic energy, vorticity 
dynamics, and stability considerations. These 


different considerations oust be consistent with 
each other and, therefore, only provide different 
viewpoints from which the flow can be studied. 

The study of the flow in terms of the vorticity 
dynamics explored the two actual mechanisms that do 
notact in constant density flows: the barocl ini c 

torque and thermal expansion. The action of these 
mechanisms was shown to result in more diffuse and 
weaker vortices when heat release accompanies the 
chemical reaction. At the largest scales (which 
dominate the dynamics of the flow and account for 
most of the turbulent transport), the altered 
vorticity distribution resulted in slower rollup of 
the most unstable modes, giving lever growth rates 
and less entrainment of unmixed fluid. This was 
indicated by the lower energies and growth rates of 
the unstable modes computed in the simulations and 
confirmed by a stability analysis of a simplified 
model problem similar to the mixing layer flow 
simulated here. The appearance of "humps" in the 
velocity profile were shown to be the result of 
vorticity generation in the outer regions of the 
vortices by barocl inic torques. 

The turbulent kinetic energy and the turbulent 
stresses are closely related, since the turbulent 
stresses are a direct indication of the kinetic 
energy transfer from the mean flow to or from the 
turbulence. Lower values of both the turbulent 
kinetic energy and the turbulent stresses were 
observed in the heat release runs. This can be 
related directly to the vorticity dynamics by 
realizing that in turbulent flows, the largest 
eddies (vortices) are responsible for most of the 
transport of momentum and scalar variables. The 
weaker large-scale vortices that result when heat 
release occurs transport less momentum, which is 
exactly what the lower values of the turbulent 
shear stresses indicate. The transport rates of 
the chemical reactants will also result in less 
product formation, as observed in the simulations 
and similar experiments. 
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TABLE 1 

THREE-DIMENSIONAL SIMULATION PARAMETERS 


Run No. 

Domain 

Re 

Da 

A i.o 

A l/2,0 

*3D 

Ce 

A 

1 

4 ft 

500 

2 

0.0 

0.0 

0.00003 

0 

2 

2 

Att 

500 

2 

0.0 

0.0 

0.00003 

5 

2 

3 

4 71 

500 

2 

0.01 

0.01 

0.00003 

0 

2 

4 

Att 

500 

2 

0.01 

0.01 

0.00003 

5 

2 


Al,0 A 1/2,0 «re the amplitude of most unstable node and its 

subharmonic as determined from a linear stability analysis of the 
hyperbolic tangent velocity profile .5° 
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Figure 5 Spanwiae conponent of vorticity. Two- 
dinentional * inula t ion, no heat, runl. 








Figure 11 Vorticity thickness, three-dimensional 
simulation, run3 (no heat), runA (heat 
release) . 


Figure 12 Velocity half width, three-dimensional 
simulation, run3 (no heat), runA (heat 
release). 
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Figure 19 Turbulent transport term in kinetic Figure 20 Velocity-pressure gradient correlation, 

energy equation, runl (no heat), run2 runl (no heat), run2 (heat release), t»48. 

(heat release), t*48. 



22a t-48. 



Figure 22 Spanwise component of vorticity. 

Three-dimensional simulation, no heat, 
runl . 


22b t-72 
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23a t-48. 


Contour From 

F 0.70 to 0.2 
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23b t-72. 


Figure 23 Spanwise component of vorticity. Three- 
dimensional simulation, heat release, 
run2. 



Figure 24 Instantaneous value of expansion term in 
vorticity equation, runl, x“0. 
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25a t-48. 25b t-72. 

Figure 25 Instantaneous value of baroclinic 

torque. Three-dimensional simulation, 
run2. 



26a t“12. 


26b t-24. 



26c t-48. 


26d t-72. 


Figure 26 Average value of baroclinic torque and 

expansion terms plotted across the mixing 
layer, run2. 
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Contour From 
-0.30 to 0.50 


Contour From 
-0.80 to 1.40 


Figure 27 Streann.’ise component of vorticity, run2 
(heat release) streatrvise location x“tt, 
t‘72. 


Figure 28 Streamwise component of vorticity, runl 
(no heat release), atreanvise location 
x"H, t«72. 
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Figure 29 Configuration for atability analysis. 
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Figure 30 Stability anaysi6 result* 
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Figure 31 Energy contained in various vsvnumbers, Figure 32 Energy contained in various wavenumbers, 

runl (no heat). run2 (heat). 
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Figure 33 Comparison of energy in fundamental node. Figure 34 Comparison of energy in subharmonic, runs 
runs 1 and 2. 1 and 2. 



as ss sc ts 

HUE 


Figure 35 Comparison of energy in three- 
dimensional modes, runs 1 and 2. 




Figure 36 Energy contained in various wavenumbers , 
run3 (no heat). 


Figure 37 Energy contained in various wavenumbers 
run4 (heat). 
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